简介:一套基于MATLAB的层次分析法(AHP)实用工具,聚焦判断矩阵权重计算与一致性检验,适合需要处理多准则决策问题的研究人员、工程师和学生,可应用于项目评估、方案优选等场景。压缩包共8个文件,以5个m脚本为核心,分别实现算术平均值法、特征值法以及CR一致性校验;同时配套1个xlsx例题数据、1个txt成对比较矩阵和1个eddx层次结构图,整套资源仅31KB,轻量易用。已有589人学习。代码允许用户直接输入自定义判断矩阵,自动完成权重求解、最大特征值λmax提取和一致性比率CR计算,并给出是否通过一致性检验的判断,帮助深入理解AHP中RI、CR等关键概念。配套的例题数据与说明文件便于对照验证,可显著简化多准则决策分析流程,适合作为课程设计、论文研究或实际项目中的参考工具。
1. 层次分析法判断矩阵权重计算的MATLAB实现
手算层次分析法权重时,最典型的翻车点不是比较矩阵填错,而是求最大特征值λmax。3阶矩阵勉强能用特征方程手算,5阶以上几乎没有手算的可能性,一旦λmax算错,一致性比率CR跟着错,决策结论就全站不住脚。这套MATLAB代码和例题数据的价值在于:把判断矩阵的读入、三种权重计算、一致性检验、方案得分排名串成了一条可以直接跑通的流程。适合正在做毕业论文的工科学生、需要做供应商评估或技术选型的工程师,以及想把AHP从Excel公式里解放出来的人。输入一个成对比较矩阵,输出权重向量、CR值和方案得分,整个过程是黑箱可验证的,而不是拿别人的结论当结论。
2. 判断矩阵构造与三种权重计算方法
2.1 判断矩阵的1-9标度与填写规则
判断矩阵是AHP的输入,它表达的是“同一层次里,元素i相对于元素j的重要程度”。标度规则只有一套:1代表同等重要,3表示前者稍微重要,5表示明显重要,7表示强烈重要,9表示极端重要,2、4、6、8用于中间状态。倒数则表达反方向比较,a_ji = 1/a_ij。对角线上的元素永远是1,因为自己和自己同等重要。
实际填矩阵的时候,最容易犯的错误是漏填左下角。因为只填上三角再取倒数补全下三角,这种操作在Excel里会留下公式依赖,一旦复制顺序出错,矩阵就变成非互反矩阵,直接导致权重计算结果不可用。我一般建议在MATLAB里用triu函数构建上三角,再通过转置倒数合成完整矩阵,而不是手工输入全部n²个元素。下面的代码演示了这种安全的构造方式,确保矩阵满足互反性。
% 构造3阶判断矩阵,只填写上三角部分 n = 3; A = zeros(n); % 上三角元素:1行2列=3,1行3列=5,2行3列=2 A(1,2) = 3; A(1,3) = 5; A(2,3) = 2; % 对角线为1 A(1,1) = 1; A(2,2) = 1; A(3,3) = 1; % 通过转置倒数补全下三角 A = A + (A' .* (A == 0));这段代码的核心逻辑是A'把上三角对称过来,. * (A == 0)保证只有原来为零的位置被填充,不会覆盖已经填好的上三角值。补全后建议用issymmetric(A .* 0 + 1)检查结构,再用any(any(A .* A' ~= 1))验证互反性。参数n代表比较元素个数,3阶矩阵只做3次比较就能确定全部6个非对角线元素,这是AHP对决策者认知负担的一种妥协。
2.2 算术平均法求权重:公式与MATLAB函数
算术平均法是最直观的一种权重算法。它的思路是:把判断矩阵的每一列归一化,得到新的列向量,每个列的归一化向量都代表一种“伪权重”,再对行求平均,得到最终权重。公式为:
```math w_i = \frac{1}{n} \sum_{j=1}^{n} \frac{a_{ij}}{\sum_{k=1}^{n} a_{kj}}
用MATLAB实现时,不需要循环,直接用矩阵除法。`sum(A,1)`返回每列的和,`A ./ sum(A,1)`利用广播机制完成列归一化,`mean(result,2)`对行求平均。下面这个函数`ahp_arithmetic`就是资源中`ccfx_base.m`里核心算法的简化版本,一行计算就完成了整个流程。 ```matlab function w = ahp_arithmetic(A) % 算术平均法求判断矩阵权重 % 输入A为n阶正互反矩阵 % 输出w为n维权重向量,元素和为1 n = size(A, 1); % 列归一化:每个元素除以所在列的和 normA = A ./ sum(A, 1); % 对归一化后的行取平均 w = mean(normA, 2); % 保证权重和为1(浮点误差修正) w = w / sum(w); end这里sum(A, 1)返回1xn的行向量,A是nxn矩阵,两者相除时MATLAB隐式扩展,等价于每一列除以对应的列和。mean(normA, 2)沿行方向求均值,得到n维列向量。最后一步w / sum(w)是预防浮点累计误差,因为列归一化和行平均理论上已经保证和为1,但浮点运算可能导致第n位小数偏移。这个方法计算速度最快,但它的理论依据是各列地位平等,实际在比较矩阵一致性较差时,算术平均法的结果会偏离特征向量法,适合快速预览权重分布。
2.3 几何平均法与特征值法的代码实现
几何平均法(又称乘积方根法)先把每一行元素相乘,再开n次方,最后归一化。它的数学表达式是 w_i = (∏ a_ij)^(1/n) / Σ (∏ a_kj)^(1/n)。用MATLAB的prod函数可以无循环实现。特征值法则是取判断矩阵最大特征值对应的归一化特征向量作为权重,这是Saaty原始论文里规定的方法,也是多数决策软件默认采用的方法。特征值法能同时给出λmax,后续一致性检验需要这个值。
function [w, lambda_max] = ahp_eig(A) % 特征值法求判断矩阵权重,同时返回最大特征值 % 输入A为n阶正互反矩阵 % 输出w为权重向量,lambda_max为最大特征值 [V, D] = eig(A); % 提取对角线上的特征值 lambda = diag(D); % 找到最大特征值的索引 [lambda_max, idx] = max(lambda); % 取对应的特征向量 v = V(:, idx); % 归一化特征向量作为权重 w = abs(v) / sum(abs(v)); endeig(A)返回特征向量矩阵V和特征值对角阵D,diag(D)把对角元素抽成向量。max返回的第二个值idx是最大特征值的位置,对应特征向量就是主特征向量。这里取了绝对值是因为矩阵可能存在负特征向量分量,但正互反矩阵的Perron-Frobenius定理保证主特征向量全为正,abs只是防御性写法。注意lambda_max在输出里是一个标量,可以直接传递给一致性检验函数。特征值法的精度最高,也是资源中ccfx.m采用的方法,但计算复杂度为O(n³),对10阶以上矩阵才开始有感知差异,决策场景完全可接受。
2.4 三种方法的对比与适用场景
| 方法 | 核心计算 | 抗一致性偏差能力 | 计算速度 | 适用场景 |
|---|---|---|---|---|
| 算术平均法 | 列归一化求行均值 | 较弱 | 最快 | 快速预览、上课演示 |
| 几何平均法 | 行乘积开n次方 | 中等 | 快 | 矩阵阶数较高时稳健 |
| 特征值法 | 特征分解 | 最强 | 较慢 | 正式决策、学术论文 |
三种方法计算同一个一致性良好的矩阵,权重误差通常在0.01以下。当矩阵一致性较差时,特征值法会放大不一致性带来的影响,这反而是优势——它把问题暴露出来。几何平均法的数学性质介于两者之间,而且永远返回正权重,不会因为数值问题出现零权重。如果只把AHP当排序工具而不做严格学术报告,算术平均法足够;如果要用CR值说服评审,必须用特征值法,因为一致性检验的理论基础就建立在最大特征值上。资源里的ccfx_Learn.m同时包含三种方法,可以用同一组数据对比输出,这种对比本身就是对算法正确性的验证。
3. 一致性检验:CI、RI、CR的完整计算
3.1 一致性指标CI与平均随机一致性RI
判断矩阵允许不一致,但必须控制在范围内。一致性指标的定义是 CI = (λmax - n) / (n - 1),这里的 n 是矩阵阶数。当 λmax = n 时,CI = 0,代表矩阵完全一致;CI 越大,不一致程度越高。但CI不能直接作为判断标准,因为随机矩阵也会产生随机的CI,所以需要引入平均随机一致性指标RI。
RI 的值来自Saaty通过500次随机模拟实验得出的统计表,它只与矩阵阶数 n 有关。1阶和2阶矩阵永远完全一致,从3阶开始需要检验。下面的表是标准RI数值,不同文献略有差异,但主流论文使用这个版本:
| n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| RI | 0 | 0 | 0.58 | 0.90 | 1.12 | 1.24 | 1.32 | 1.41 | 1.45 | 1.49 |
注意RI值针对的是正互反矩阵,不是任意实数矩阵。如果你的判断矩阵是残缺的(比如允许空值),RI需要另行计算,不能用这个表。在MATLAB里我习惯把RI存成常量向量,按阶数索引直接取值,而不是用switch-case,这样代码更紧凑。
3.2 一致性比率CR及其MATLAB实现
一致性比率 CR = CI / RI,当 CR < 0.1 时,认为判断矩阵的一致性可接受;否则需要修正。这个0.1阈值是经验值,Saaty当初提出时没有严格统计推导,但业界一直沿用。下面是完整的一致性检验函数,结合了3.1中的RI表。
function [CR, CI, lambda_max] = ahp_consistency(A) % 判断矩阵一致性检验 % 输入A为n阶正互反矩阵 % 输出CR为一致性比率,CI为一致性指标,lambda_max为最大特征值 n = size(A, 1); % 特征值法计算权重和最大特征值 [~, lambda_max] = ahp_eig(A); % 一致性指标 CI = (lambda_max - n) / (n - 1); % 平均随机一致性指标表 RI_list = [0, 0, 0.58, 0.90, 1.12, 1.24, 1.32, 1.41, 1.45, 1.49]; % 1阶和2阶矩阵无需检验 if n <= 2 CR = 0; else RI = RI_list(n); CR = CI / RI; end % 输出是否通过检验 if CR < 0.1 fprintf('一致性检验通过,CR = %.4f\n', CR); else fprintf('一致性检验失败,CR = %.4f,请修正判断矩阵\n', CR); end endahp_eig返回两个值,用~丢弃第一个权重输出。RI_list(n)用矩阵阶数作为索引,如果n超过10,这个函数会越界,实际中超过10阶的判断矩阵很少见,但严谨做法是在函数开头加assert(n <= length(RI_list))。输出使用fprintf是为了方便在命令行看到结果,不是必须的。注意这里CR计算用的是浮点除法,CI有可能因为λmax的数值误差出现极小负数,比如-0.0002,这种情况应该手动置为0,否则CR为负值会被误判为完全一致。可以在CI计算后加一行CI = max(CI, 0);。
3.3 检验失败时如何定位问题元素
CR超过0.1时,多数人第一反应是全部重填矩阵,这既浪费时间又未必有用。更好的做法是找到最“不一致”的那个比较值。我常用的手段是:先用当前权重向量w重建一个理想一致矩阵B,B_ij = w_i / w_j,然后计算原始矩阵A与B的差异矩阵D = A - B,差异绝对值最大的位置就是问题元素。
function [i, j, max_diff] = find_inconsistency(A, w) % 定位判断矩阵中最不一致的元素位置 % 输入A为判断矩阵,w为权重向量 % 输出i,j为问题元素下标,max_diff为最大绝对差异 B = w * (1 ./ w)'; % 理想一致矩阵 D = abs(A - B); % 只考察上三角,避免重复报告 D = triu(D, 1); [max_diff, lin_idx] = max(D(:)); [i, j] = ind2sub(size(A), lin_idx); fprintf('问题元素 A(%d,%d)\n', i, j); endw * (1 ./ w)'是外积运算,产生n阶矩阵,每行是wi除以所有wj,正好符合一致矩阵的定义。triu(D,1)取严格上三角是因为下三角是上三角的倒数,差异模式相同,避免重复报告。ind2sub把线性索引转换成行和列。执行这个函数后,通常会发现最大差异集中在某几个元素上,把对应标度往1的方向调整(比如5改成3,3改成2),重新计算CR就能通过。注意不要同时修改多个元素,否则很容易从一种不一致跳到另一种不一致。
4. 例题数据实战:从成对比较矩阵到得分排名
4.1 读取“成对比较矩阵.txt”中的判断矩阵
资源包里的成对比较矩阵.txt是纯文本格式,每行对应矩阵一行,元素用空格分隔。用MATLAB的readmatrix可以一次性读入,但要注意文件里不能有中文表头或注释行。如果文件里有分隔符不统一的情况,我一般会用readmatrix('成对比较矩阵.txt', 'Delimiter', ' ', 'NumHeaderLines', 0)显式指定空格作为分隔符。读取后先用size确认矩阵阶数,再调用前面的ahp_consistency做一次快速体检,避免后续计算基于一个不一致的输入。
% 读取判断矩阵文件 A = readmatrix('成对比较矩阵.txt'); % 显示矩阵以确认格式正确 disp(A); % 进行一致性预检,提前发现数据问题 [CR, ~, ~] = ahp_consistency(A); if CR >= 0.1 warning('判断矩阵一致性不达标,检查数据后再继续'); endreadmatrix是MATLAB R2019a引入的函数,老版本里对应的是dlmread。如果文件最后一列带换行符异常,可能出现读取结果维度不对,这时候用load('成对比较矩阵.txt')也能读,但要求文件是纯数字矩阵。建议读取后立刻加一个断言assert(size(A,1) == size(A,2))。如果判断矩阵是8阶以上,建议检查A中是否有负数或0,因为正互反矩阵只允许正数。
4.2 计算方案层得分并写入Excel
在“根据权重矩阵计算得分.xlsx”里,通常保存着一组方案评分数据:行是备选方案,列是评价指标,每个单元格是某个方案在某个指标下的得分。AHP最终的决策得分是权重向量与评分矩阵的乘积。这里的关键是确保评分矩阵的列数等于权重向量的维度,否则矩阵乘法会报维度不匹配。
% 读取方案评分矩阵score_data score = readmatrix('根据权重矩阵计算得分.xlsx'); % 用特征值法计算指标权重 w = ahp_eig(A); % 这里w是n维列向量 % 计算综合得分 = score * w final_score = score * w; % 展示结果 results = table((1:size(score,1))', final_score, ... 'VariableNames', {'方案编号', '综合得分'}); disp(results); % 写入结果到新Excel writetable(results, '得分结果.xlsx');score * w是标准的矩阵乘法,每个方案的得分等于该方案在各指标得分加权求和。如果score的列数不等于w的长度,MATLAB会直接报错“矩阵维度必须一致”。这里要注意,score矩阵的每一列对应一个判断矩阵中的元素,而不是整列代表一个指标的全部方案得分。在决策实战里,评分数据往往有量纲差异,建议先对评分矩阵做min-max归一化再乘权重,否则量纲大的指标会主导得分。资源里的ccfx.m还包含一个把得分排序并输出排名的片段,常见做法是用sort(final_score, 'descend')得到降序索引,再用排名=1:length(final_score)构造排名列。
4.3 用亿图eddx核对层次结构
层次分析法.eddx是亿图图示(Edraw Max)绘制的层次结构图文件,不参与数值计算,但它承载了决策模型的目标层、准则层和方案层关系。当判断矩阵的阶数超过5时,单看矩阵很难发现层次结构错误,比如把准则层元素和方案层元素放进了同一个判断矩阵。我一般会把eddx打开后,对照层次图检查矩阵的行列顺序。如果发现顺序不一致,需要修改矩阵行列而非权重代码,因为代码只是按矩阵行号输出权重,不会自动识别元素名称。
一种保证顺序一致的方法是:在MATLAB里把指标名称定义成字符串向量,并在读取矩阵时同时读取列名。如果eddx里的层次顺序是“价格、性能、售后服务”,那么判断矩阵的行列顺序也必须是这个顺序。这里有一个常见错误:计算得分时,评分矩阵的列顺序与判断矩阵的行顺序不一致,导致最终权重乘错列。为了避免这个问题,建议在写入得分结果时,把指标名称一起写到表头里,而不是只用“指标1、指标2”。writetable支持自定义列名,用'VariableNames'参数即可。
5. 提高判断矩阵权重计算效率的批量处理技巧
5.1 批量计算多个判断矩阵的无循环写法
决策场景经常有多个专家各自打分,产生一批判断矩阵。逐个调用ahp_eig没有错,但更简洁的做法是把所有矩阵堆叠成一个三维数组,用arrayfun配合匿名函数批量计算。假设有k个专家的判断矩阵,每个都是n阶方阵,把它们放在A_all的第三维里。
% A_all是n x n x k的三维数组 n = size(A_all, 1); k = size(A_all, 3); % 批量计算权重 W = zeros(n, k); Lambda = zeros(k, 1); for t = 1:k [W(:,t), Lambda(t)] = ahp_eig(A_all(:,:,t)); end % 输出全部权重,每列对应一个专家 disp(W);循环在这里比arrayfun更清晰,因为ahp_eig有两个返回值,arrayfun需要额外把两个输出合并成一个元胞数组。参数t是专家序号,A_all(:,:,t)取出第t个矩阵。批量计算的瓶颈不在循环次数,而在每次eig调用的矩阵阶数。如果同时计算20个5阶矩阵,总耗时不到0.05秒,完全不需要优化。
5.2 用CR的贡献量辅助矩阵修正
修正不一致矩阵时,单纯看差异矩阵的绝对值可能不够,因为某些位置的差异对CR贡献大,另一些则无关紧要。更精确的做法是计算每个元素变化对CR的影响量。这个方法思路是:把某元素a_ij替换成1(即认为i和j同等重要),重新计算CR,对比CR的变化幅度。变化大的位置就是修正优先级高的位置。
function sensitivity = cr_sensitivity(A) % 计算每个位置元素变为1后对CR的影响量 % 输出sensitivity为n阶矩阵,仅上三角有值 n = size(A, 1); base_CR = ahp_consistency(A); sensitivity = zeros(n); for i = 1:n-1 for j = i+1:n temp = A; temp(i,j) = 1; temp(j,i) = 1; [CR, ~, ~] = ahp_consistency(temp); sensitivity(i,j) = abs(CR - base_CR); end end endbase_CR是原始矩阵的CR值,temp是替换后的矩阵,abs(CR - base_CR)量化了该元素对整体一致性的影响。如果某个位置的sensitivity很大,说明它是“关键不一致源”。这里的嵌套循环在n=5时会执行10次,时间可忽略。注意ahp_consistency会打印CR值,批量调用时命令行会被刷屏,可在函数内部加一个silent参数抑制输出。
5.3 把权重结果整理成可复现的表格
权重算完后,最容易被评审挑毛病的是缺少可复现性。我会把原始判断矩阵、权重、CR一起写入同一个工作簿的不同sheet,方便事后核对。用writetable分别写矩阵和权重,再用writecell写标题行。
% 汇总输出 T_weight = array2table(w', 'VariableNames', {'权重'}); T_matrix = array2table(A); % 写入Excel工作簿 writetable(T_matrix, 'AHP报告.xlsx', 'Sheet', '判断矩阵'); writetable(T_weight, 'AHP报告.xlsx', 'Sheet', '权重'); % 记录CR值到单元格 CR_sheet = cell(1,1); CR_sheet{1,1} = 'CR值'; CR_sheet{2,1} = CR; writecell(CR_sheet, 'AHP报告.xlsx', 'Sheet', '一致性检验', 'Range', 'A1');array2table把矩阵转成表格,VariableNames是自动生成的列名,默认用“A1、A2”这类名称,建议改成实际指标名。writecell的Range参数指定写入区域,'A1'表示从左上角开始。这样做的好处是,Excel里的数据格式固定,任何人打开都能复现计算路径。最后提醒一个细节:writetable在覆盖同名文件时会直接覆盖不提示,如果循环里多次写入,前面的数据会丢失,建议文件名带上时间戳,用datestr(now, 'yyyymmdd_HHMM')生成后缀。
本文还有配套的精品资源,点击获取