用Matlab實現(xiàn)AHP的算法

1、1MATLAB 的基本內(nèi)容MATLAB（MATrix LABoratory，矩陣實驗室的縮寫） 是一種特殊用途的計算機程序優(yōu)化執(zhí)行工程和科學計算。 它開始為旨在執(zhí)行矩陣數(shù)學程式的生活，但多年來它已發(fā)展成為一個靈活的計算系統(tǒng)基本上能夠解決任何技術(shù)問題。MATLAB具有編程語言的基本特征，使用MATLAB也可以使用像BASIC、FORTRAN、C 等傳統(tǒng)編程語言一樣，進行程序設計，而且簡單易學、編程效率高。正因為 MATLAB的強大的功能， 使得它在許多領(lǐng)域得到廣泛應用。在科研與工程應用領(lǐng)域， MATLAB已被廣泛地用于科學研究和解決各種具體的實際問題。許多科技工作者選用 MATLAB做為計算工具

2、， 避免了繁瑣的底層編程， 從而可以把主要精力和時間花在科學研究和解決實際問題是上，提高了工作效率。MATLAB矩陣矩陣是 MATLAB的基本處理對象，因此根據(jù)本文所需，簡單介紹所涉及MATLAB矩陣內(nèi)容。MATLAB矩陣的建立1、直接輸入法最簡單的建立矩陣的方法是從鍵盤直接輸入矩陣的元素。例如：A=1 2 3;4 5 6;7 8 9A =123456789也可以用回車鍵代替分號，按下列方式輸入：A=123456789 2、利用 M 文件建立矩陣比較大且復雜的矩陣，可以為它專門建立一個M 文件，如同下例。利用 M 文件建立矩陣。啟動有關(guān)編輯程序或MATLAB文本編輯器，并輸入待建矩陣：MYMA

3、T= 111 ,112,113,114,115,116,117,118,119;211,212,213,214,215,216,217,218,219;把輸入的內(nèi)容以純文本方式存盤（設文件名為）。在 MATLAB 命令窗口中輸入mymatrix，即運行該M 文件，就會自動建立一個名為MYMAT 的矩陣，可供以后使用。矩陣的特征值與特征向量特征值和特征向量在科學研究和工程計算中都有非常廣泛地應用。在MATLAB中，計算矩陣 A 的特征值和特征向量的函數(shù)是 eig(A),常用的調(diào)用格式有3 種14：E = eig( A ) ：求矩陣A 的全部特征值，構(gòu)成向量E。V,D=eig(A)：求矩陣A 的全

4、部特征值，構(gòu)成對角矩陣D，并求A 得特征向量構(gòu)成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 的各列構(gòu)成的向量。MATLAB的 M 文件用 MATLAB語言編寫的程序稱為 M 文件。 M 文件是由若干 MATLAB命令組成在一起構(gòu)成的，它可以完成某些操作，也可以實現(xiàn)某種算法。M 文件可以根據(jù)調(diào)用方式的不同分為兩類：命令文件（Scrip

5、t File）和函數(shù)文件（FunctionFile）。它們的擴展名均為 .m。函數(shù)文件由 function 語句引導，其基本結(jié)構(gòu)為：function輸出形參表 =函數(shù)名（輸入形參表）注釋說明部分函數(shù)體語句我們通過舉例說明如下：例 2-2 分別建立命令文件和函數(shù)文件， 將求矩陣的一致性指標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 建立函數(shù)文件。function c=C

6、I(max,n)c=(max-n)/(n-1)然后在 MATLAB的命令窗口調(diào)用該函數(shù)文件。max=input(please input max:);n=input(please input n:);c=CI(max,n)2基于 MATLAB的 AHP實現(xiàn)21 AHP的 MATLAB的計算流程框圖根據(jù)層次分析法的一般步驟我們得到在圖，如圖 2 所示 16開 始輸入準則層層數(shù)輸入準則層對目標層aij判斷矩陣元素判斷矩陣歸一化計算R0.1輸 入 準 則 層 第 j 個 準 則包 含 的 方 案 數(shù) 量 n輸 入 方 案 層 第 j 個 準 則層 的 判 斷 矩 陣 元 素 bMATLAB 工具上實

7、現(xiàn)的計算程序流程框結(jié) 束CR0.1總CR計算Fa權(quán)值的計算計算權(quán)值矩陣輸 入 方 案 層 第 m 個 準則層的相關(guān)矩陣元素判斷矩陣歸一化計算C R 0(2) aij a ji = l(3) aii =l按照事物邏輯要求，該矩陣還應具備一致性，即滿足：aij a jk = aik前面已經(jīng)給出由于客觀事物的復雜性與決策者的認識的多樣性，實際問題的成對比較矩陣不可能做到嚴格上的一致性，因而，借助平均隨機一致性指標RI 來相對判定其一致性程度。其中表1-4 是已經(jīng)計算好的 115 階矩陣的 RI 值表，但未給出其實現(xiàn)過程，且各文獻的RI 值表不完全相同。究其原因除沒有太大必要介紹外，真正去實現(xiàn)它卻有如

8、下三個難度 ：(1)隨機兩兩判斷矩陣中的元素要求是19 和它們對應的倒數(shù)共17 個整數(shù)與小數(shù)的均勻分布很難處理。(2)一般高級編程語言實現(xiàn)成對比較矩陣及相關(guān)計算，非常復雜，且占用內(nèi)存巨大，耗時多。隨機種子源不能控制。本文使用數(shù)學軟件包MATLAB對其進行計算。 其設計解決思路為： 先用軟件包中隨機函數(shù)產(chǎn)生數(shù)117 的均勻分布的 n 階矩陣，然后在軟件包中采用不同技巧將它轉(zhuǎn)化為成對比較矩陣，最后用循環(huán)語句計算出RI 值。結(jié)果如下表：表 3-1 計算的 RI 值階數(shù)123456789RI00以 n=4 為例，過程詳見如下程序清單， 其中隨機成對比較矩陣的實現(xiàn)見相應注釋部分。 15MATLAB的程序

9、 M 文件：function rin=4; ri=0; m=100;%計算RI 值的命令文件rand(seed,21)for i=1:m%控制隨機發(fā)生器a=ceil (17*rand(n);%產(chǎn)生n 階l17 的隨機陣a(find(a=8)=;%消除0 為分母b=1./(a-8);%產(chǎn)生一個輔助陣a(find(a9)=b(find(a9);%借助b，將917 分別轉(zhuǎn)化為1129a(find(a=)=8;e=eye(n);%產(chǎn)生一個4 階單位陣c=1./a;%將a 中每個元素換成相應倒數(shù)c=c;%將c 轉(zhuǎn)置c=tril(c,-1);%抽取c 的下三角(不含主對角線)a=triu(a,1);%抽取

10、a 的上三角 (不含主對角線)a=a+c+e;%實現(xiàn)隨機成對比較陣ak=size(a,1);%計算a 的行維數(shù)ri=ri+(max(abs(eig(a)-k)/(k-1);%計算100 次RI 值endri/m%計算平均RI 值23 AHP各環(huán)節(jié)的MATLAB實現(xiàn)1351135以目標矩陣 A= 1 313213,準則層矩陣為 P1=, P2= 13，1 51 311 211 311 512P3=為例，運用 MATLAB進行數(shù)據(jù)處理。特征向量及其歸一化的MATLAB實現(xiàn)MATLAB中求矩陣特征值和特征向量的函數(shù)是eig，其調(diào)用的格式為 V,D=eig(A),其中， V 為特征向量矩陣， D 為特

11、征值矩陣。層次分析法中需要求得是最大特征值及對應的歸一化特征向量，而且考慮到eig 函數(shù)在求得的特征值中可能會存在復數(shù)。因此，運用直接輸入程序代碼會產(chǎn)生一定的誤差。在此需要對求得的V、D 進行適當選擇，定義一個 M-file 來實現(xiàn)。functionmaxeigval,w=maxeigvalvec(A)%求最大特征值及對應的歸一化特征向量%A 為判斷矩陣eigvec,eigval=eig(A);eigval=diag(eigval);%特征向量eigvalmag=imag(eigval);realind=find(eigvalmag=disp(input( 矩陣沒通過一致性檢驗，請重新調(diào)整判斷

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);運行結(jié)果如下：矩陣通過一致性檢驗矩陣通過一致性檢驗矩陣通過一致性檢驗矩陣通過一致性檢驗一致性檢驗及總排序的MATLAB實現(xiàn)通過層次單排序（權(quán)重）計算后，進行層次合成計算，在此本文定義函數(shù)計算層次總排序的權(quán)重并進

13、行一致性檢驗。function tw=tolsortvec(utw,dw,CIC,RIC)求層次總排序權(quán)重并進行一致性檢驗utw 為上一層因素的總排序權(quán)重行向量dw 為下一層因素相對于上一層各因素的層次單排序權(quán)重矩陣CIC為一致性指標列向量RIC為隨機一致性指標列向量tw=dw*utwCR=utw*CIC/(utw*RIC);if CR=disp(input( 層次總排序沒通過一致性檢驗，請重新調(diào)整判斷矩陣);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);運行結(jié)果如下：tw=層次總排序通過一致性檢驗其中 tw 是層次總排序結(jié)果。因此，根據(jù)數(shù)據(jù)建立如下的層次總排序表。表 3-3層次總排序表（權(quán)重）AA層次P的注總排序結(jié)果PP10000P0020000P30000注：按概率乘法，P 層次總排序指標的權(quán)重值為N P 層次指標的權(quán)重與相應上一層次指標A N 層權(quán)重的積，且總排序權(quán)重值的和為1。選擇最優(yōu)排序計算出層次總排序后， 為了使決策者能迅速得出結(jié)果， 本文對層次總排序進行最優(yōu)排序。

15、運用 MATLAB鍵入如下指令：n=length(tw);for i=1:nt=max(tw);b(i)=t;m n=find(a=t);tw(n)=;endb運行結(jié)果如下：b=利用 MATLAB大大縮短了計算復雜矩陣的時間，為決策者節(jié)省了寶貴的時間， 從而有更多的精力投入其他事務。3基于 MATLAB的 AHP應用31挑選合適工作問題某畢業(yè)生選擇工作，經(jīng)雙方懇談，假設已有三個單位C1,C2,C3 表示愿意錄用他。該生對三個單位進行了解后，選取了一些中間指標進行考察，例如單位的研究課題，發(fā)展前途，待遇，同事情況，地理位置，單位名氣等。根據(jù)層次分析法，試求該生工作優(yōu)先排序（給出權(quán)值、計算程序）

16、，并給出最終選擇決策?，F(xiàn)以 A、B、C 表示選擇工作的三個層次，建立如下結(jié)構(gòu)模型：圖 3 選擇單位層次結(jié)構(gòu)圖根據(jù)成對比較法，得到相應判斷矩陣如下表：表 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現(xiàn)在在 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( 準則層對目標層權(quán)向量); disp(WA);disp( 方案層對準則層權(quán)向量); disp(WK);E=WK*WAdisp( 方案層組合權(quán)向量);

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運行結(jié)果如下：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通過一致性檢驗!準則層對目標層權(quán)向量方案層對準則層權(quán)向量E =方案層組合權(quán)向量組合一致性比率CR通過一致性檢驗!CHOICE =12）、M 文件法：根據(jù)第三章節(jié)的層次分析法各環(huán)節(jié)的MATLAB實現(xiàn)定義如下三個函數(shù)：（1） 定義來實現(xiàn)最大特征根及對應的歸一化特征向量。functionmaxeigval,w=maxeigvalvec(A)%求最大特征值及對應的歸一化特征向量%A 為判斷矩陣eigvec,eigval=eig(A);eigval=diag(eigval);%特

24、征向量eigvalmag=imag(eigval);realind=find(eigvalmag=disp(input( 矩陣沒通過一致性檢驗，請重新調(diào)整判斷矩陣)elsedisp(input( 矩陣通過一致性檢驗 );end（3）定義函數(shù)計算層次總排序的權(quán)重并進行一致性檢驗。function tw=tolsortvec(utw,dw,CIC,RIC)求層次總排序權(quán)重并進行一致性檢驗utw 為上一層因素的總排序權(quán)重行向量dw 為下一層因素相對于上一層各因素的層次單排序權(quán)重矩陣CIC為一致性指標列向量RIC為隨機一致性指標列向量tw=dw*utw;CR=utw*CIC/(utw*RIC);if

25、CR=disp(input( 層次總排序沒通過一致性檢驗，請重新調(diào)整判斷矩陣);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