1. 项目缘起为什么我们需要一份“活”的数学建模代码笔记如果你也参加过数学建模竞赛或者正在学习相关课程大概率会和我有一样的感受网上的代码资源看似丰富实则“坑”多“雷”也多。你兴冲冲地下载了一个号称“层次分析法AHP完整代码”的压缩包解压后却发现要么是几行语焉不详的脚本注释比代码还少要么是依赖某个早已过时的第三方库运行起来满屏报错更常见的是代码逻辑僵硬只能处理教科书上的理想案例稍微变一下数据格式或者评价指标就彻底歇菜。这种经历我称之为“代码废墟考古”耗费大量时间在环境配置和调试上真正用于理解模型、优化方案的时间所剩无几。我最初整理这份“清风数学建模代码笔记”的动机正是源于无数次这样的挫败。我不想要那种冷冰冰的、只能“看”不能“用”的代码陈列。我想要一份“活”的笔记——这里的“活”意味着每一行代码都经过实际数据的检验每一个函数都配有清晰的输入输出说明和边界条件提示每一个算法实现都附带了“为什么这么写”的思考过程以及“我在这里踩过什么坑”的实战心得。这份笔记的起点正是数学建模中最经典、应用也最广泛的两个评价类模型层次分析法AHP和TOPSIS法。它们不仅是国赛、美赛的常客更是许多课程大作业、科研评价的基石。但很多教学材料止步于理论公式将矩阵计算、一致性检验、权重合成这些关键步骤留作“课后练习”让初学者望而却步。因此这份笔记的第一个目标就是彻底打通从理论到可运行代码的“最后一公里”。我将以“清风数学建模”课程的正课内容为骨架但填充进去的全是经过实战打磨的“血肉”。我们会从最基础的环境搭建讲起确保你的Python环境能无障碍运行所有代码然后我们会像拆解一台精密仪器一样把AHP和TOPSIS的每一个步骤都用代码和注释清晰地呈现出来并重点讲解那些容易出错和产生误解的细节。比如层次分析法中判断矩阵怎么录入才高效一致性检验通不过怎么办TOPSIS法里的“最优解”和“最劣解”到底怎么定义熵权法计算权重时遇到指标值全部相同导致对数无穷大的情况该如何处理这些才是实战中的真问题。这份笔记适合所有被数学建模代码困扰的朋友无论你是初次接触建模的本科生还是需要快速复现模型的研究生抑或是想寻找可靠代码模板的参赛队员。我希望它能成为你手边一份随时可以查阅、修改和信任的“代码工具箱”而不仅仅是一份阅读材料。让我们跳过那些华而不实的理论堆砌直接进入能解决实际问题的代码世界。2. 环境奠基构建一个稳定且可复现的Python建模工作流在开始敲下任何一行模型代码之前建立一个干净、稳定、易于管理的Python环境是至关重要的一步。这一步常常被忽略但却是后续所有工作能否顺利进行的基石。一个混乱的环境会导致库版本冲突、代码在自己电脑上能跑在别人电脑上就报错等一系列令人头疼的问题。我将分享一套我经过多个项目验证的、基于conda和pip的最佳实践。首先我强烈建议使用Anaconda或Miniconda来管理你的Python环境。它们可以为你创建独立的虚拟环境避免不同项目之间的依赖打架。假设你已经安装了conda我们打开终端或Anaconda Prompt创建一个专用于数学建模的新环境conda create -n math_modeling python3.9 -y这里我选择了Python 3.9这是一个在稳定性和库兼容性之间取得很好平衡的版本。当然选择3.8或3.10也可以但尽量避免使用太新或太旧的版本以免某些科学计算库尚未适配。创建完成后激活这个环境conda activate math_modeling你会看到命令行提示符前面变成了(math_modeling)这表示你已经进入了这个独立的环境。接下来我们安装最核心的几个库。我不建议一次性安装anaconda这个巨大的元包而是按需安装保持环境精简。pip install numpy1.23.5 pandas1.5.3 matplotlib3.6.3 scipy1.10.0为什么指定版本这是为了“可复现性”。今天写的代码半年后因为pandas某个函数API变了而跑不起来这种情况太常见了。锁定主要库的版本能极大避免这种问题。numpy和pandas是数据处理的基石matplotlib用于可视化scipy则提供了许多高级数学、科学计算函数在后续的矩阵运算和优化中会用到。对于数学建模还有一个非常重要的库sklearnscikit-learn。虽然它主打机器学习但其数据预处理模块如标准化、归一化和评估模块极其好用且代码经过高度优化。pip install scikit-learn1.2.2安装完成后我习惯创建一个requirements.txt文件来记录当前环境的所有依赖方便自己未来重装或分享给队友。pip freeze requirements.txt这个文件里会列出所有库及其精确版本。队友拿到你的代码和这个文件后只需要在他的环境中执行pip install -r requirements.txt就能一键复现完全相同的环境从根本上杜绝“在我这好好的到你那就出错”的尴尬。注意如果你在安装某些库时遇到速度慢或超时的问题可以临时使用国内的镜像源例如清华源pip install -i https://pypi.tuna.tsinghua.edu.cn/simple some-package。但请注意在生成requirements.txt和最终部署时应使用官方源以确保稳定性。最后我们来规划一下项目目录结构。一个清晰的结构能让你事半功倍。我通常这样组织math_modeling_notes/ ├── data/ # 存放原始数据和生成数据 │ ├── raw/ # 原始数据只读 │ └── processed/ # 清洗处理后的数据 ├── docs/ # 项目说明、参考文献 ├── notebooks/ # Jupyter Notebook用于探索性分析和演示 ├── src/ # 源代码目录 │ ├── models/ # 核心模型实现如ahp.py, topsis.py │ ├── utils/ # 工具函数如数据加载、可视化 │ └── config.py # 配置文件如路径常量 ├── outputs/ # 模型输出结果、图表 ├── requirements.txt # 依赖列表 └── README.md # 项目总览在src/models目录下我们将创建本笔记的核心代码文件。使用专业的IDE如VSCode、PyCharm打开这个项目根目录并将解释器设置为刚刚创建的math_modeling环境。至此一个坚固、专业的代码地基就打好了。接下来我们就可以放心地在上面搭建我们的模型大厦了。3. 层次分析法AHP代码实现从矩阵构建到权重计算的完整闭环层次分析法Analytic Hierarchy Process, AHP的核心思想是将复杂决策问题分解为目标、准则、方案等层次通过两两比较构造判断矩阵进而计算各元素的相对权重。理论教材会花大量篇幅讲解1-9标度法、特征向量法求权重、一致性检验公式。但在代码层面我们需要解决一系列更具体的问题如何优雅地输入判断矩阵如何高效稳定地计算特征值和特征向量一致性检验不通过时程序该如何处理或给出怎样的提示下面我将一步步拆解并给出工业级强度的代码实现。3.1 判断矩阵的输入与验证告别手动输入的低效与错误首先我们创建一个ahp.py文件。判断矩阵的输入是第一步也是容易出错的一步。我们设计一个函数它应该能接受两种输入方式一是直接传入一个二维的numpy数组二是提供一个便捷的交互方式让用户按提示输入。import numpy as np from typing import Union, Tuple def get_judgment_matrix(n: int, matrix: Union[np.ndarray, None] None) - np.ndarray: 获取或创建AHP判断矩阵。 参数: n: 矩阵维度准则或方案的数量。 matrix: 可选的预先构建好的numpy二维数组。如果提供则直接进行验证和返回。 返回: 一个符合AHP要求的numpy二维判断矩阵。 if matrix is not None: # 验证提供的矩阵 if matrix.shape ! (n, n): raise ValueError(f提供的矩阵形状{matrix.shape}与预期维度({n},{n})不符。) if not np.allclose(matrix, matrix.T): # 检查对称性 # AHP要求判断矩阵是正互反矩阵即a_ij 1 / a_ji # 这里我们检查转置后是否近似等于逆对应元素互为倒数 # 更严格的检查是对于所有i,j, matrix[i, j] * matrix[j, i] ≈ 1 for i in range(n): for j in range(i1, n): if not np.isclose(matrix[i, j] * matrix[j, i], 1.0, atol1e-9): raise ValueError(f提供的矩阵在位置({i},{j})和({j},{i})不满足互反性。) # 如果不严格对称但满足互反我们可以将其对称化取几何平均 print(警告提供的矩阵不对称将自动进行对称化处理取几何平均。) for i in range(n): for j in range(i1, n): geo_mean np.sqrt(matrix[i, j] * matrix[j, i]) matrix[i, j] geo_mean matrix[j, i] 1.0 / geo_mean return matrix.astype(float) else: # 交互式输入 print(f请根据1-9标度法两两比较以下 {n} 个元素的重要性。) print(标度含义: 1-同等重要, 3-稍微重要, 5-明显重要, 7-强烈重要, 9-极端重要, 2/4/6/8为中间值) mat np.ones((n, n)) for i in range(n): for j in range(i1, n): while True: try: value float(input(f元素 {i1} 相对于元素 {j1} 的重要性标度 (输入小数如 0.333 表示后者更重要): )) if value 0: print(重要性标度必须为正数请重新输入。) continue mat[i, j] value mat[j, i] 1.0 / value break except ValueError: print(输入无效请输入一个数字。) return mat这个函数的关键点在于鲁棒性。它允许用户直接传入矩阵方便批量处理或从文件读取也支持交互式输入。对于传入的矩阵它进行了维度和互反性的基本验证。如果矩阵不对称但满足互反性即a_ij * a_ji 1函数会给出警告并自动将其对称化这是一种常见的容错处理因为人在填写时可能因疏忽导致上下三角不完全对称。3.2 权重计算与一致性检验特征向量法的数值实现与细节处理得到判断矩阵后下一步是计算权重。教科书上会说“计算判断矩阵的最大特征值对应的特征向量并归一化得到权重”。但用代码实现时我们需要选择数值稳定的方法。numpy.linalg.eig可以计算特征值和特征向量但对于正互反矩阵其最大特征值是实数且对应的特征向量所有分量同号Perron-Frobenius定理我们可以利用这一点。def calculate_weights_and_consistency(matrix: np.ndarray) - Tuple[np.ndarray, float, float, bool]: 计算AHP判断矩阵的权重向量并进行一致性检验。 参数: matrix: AHP判断矩阵。 返回: weights: 归一化后的权重向量 (numpy array)。 lambda_max: 矩阵的最大特征值。 ci: 一致性指标 (Consistency Index)。 is_consistent: 布尔值表示是否通过一致性检验 (CR 0.1)。 n matrix.shape[0] # 计算特征值和特征向量 eigenvalues, eigenvectors np.linalg.eig(matrix) # 找到最大特征值实数部分的索引 max_eig_idx np.argmax(eigenvalues.real) lambda_max eigenvalues[max_eig_idx].real # 获取对应的特征向量取实数部分 max_eig_vec eigenvectors[:, max_eig_idx].real # 归一化特征向量得到权重 weights max_eig_vec / np.sum(max_eig_vec) # 计算一致性指标 CI (λ_max - n) / (n - 1) ci (lambda_max - n) / (n - 1) # 随机一致性指标 RI标准值这里取常用值对于n1~15 ri_dict {1: 0, 2: 0, 3: 0.52, 4: 0.89, 5: 1.12, 6: 1.26, 7: 1.36, 8: 1.41, 9: 1.46, 10: 1.49, 11: 1.52, 12: 1.54, 13: 1.56, 14: 1.58, 15: 1.59} ri ri_dict.get(n, 1.60) # 对于n15可以近似或查表这里简单处理 # 计算一致性比率 CR CI / RI cr ci / ri if ri ! 0 else float(inf) # 判断是否通过一致性检验 (通常要求 CR 0.1) is_consistent cr 0.1 return weights, lambda_max, ci, is_consistent, cr这里有几个极易踩坑的细节特征值的复数问题np.linalg.eig返回的特征值可能是复数即使对于实矩阵。由于判断矩阵是实对称的经过我们处理最大特征值应为实数。我们通过.real取实部来避免后续计算错误。特征向量的符号不确定性特征向量方向是不确定的v和-v都是对应同一特征值的特征向量。但AHP要求权重为正数。幸运的是对于正矩阵其Perron特征向量的所有分量同号。我们的代码取了实部并且后续做了归一化通常能得到正权重。如果出现负权重极少数情况说明矩阵可能不满足正互反性或者数值计算误差较大需要检查原始矩阵。RI值的选取随机一致性指标RI有标准表但不同文献略有差异。我们采用Saaty经典论文中常用的值。对于维度大于15的情况需要查阅更专业的资料或通过模拟计算得到。一致性检验的逻辑返回结果中不仅给出了布尔值is_consistent还返回了CR值本身。这是因为在实际应用中有时即使CR略大于0.1比如0.12如果决策者认为判断矩阵合理也可以接受。程序应该提供信息而不仅仅是二元判断。3.3 实战案例供应商选择决策的完整代码演示让我们用一个完整的例子来串联上述函数。假设我们要从3个供应商S1, S2, S3中选择一个评价准则有3个质量C1、价格C2、交货期C3。首先我们构建准则层的判断矩阵然后对每个准则构建方案层的判断矩阵。def ahp_supplier_selection(): AHP方法选择供应商的完整示例。 print( AHP供应商选择决策示例 ) # 1. 准则层判断矩阵 (决策者填写) # C1:质量, C2:价格, C3:交货期 # 假设决策者认为质量比价格明显重要(5)质量比交货期稍微重要(3)价格比交货期稍微不重要(1/3) criteria_matrix np.array([ [1, 5, 3], [1/5, 1, 1/3], [1/3, 3, 1] ]) print(准则层判断矩阵:) print(criteria_matrix) # 计算准则权重 criteria_weights, lambda_max_c, ci_c, consistent_c, cr_c calculate_weights_and_consistency(criteria_matrix) print(f\n准则层权重: {criteria_weights}) print(f最大特征值 λ_max: {lambda_max_c:.4f}) print(f一致性指标 CI: {ci_c:.4f}) print(f一致性比率 CR: {cr_c:.4f}) print(f一致性检验是否通过 (CR0.1): {consistent_c}) if not consistent_c: print(警告准则层判断矩阵一致性未通过请考虑调整判断值) # 2. 方案层供应商相对于每个准则的判断矩阵 # 对于准则C1质量S1比S2稍微重要(3)S1比S3强烈重要(7)S2比S3明显重要(5) matrix_c1 np.array([ [1, 3, 7], [1/3, 1, 5], [1/7, 1/5, 1] ]) # 对于准则C2价格价格越低越好所以这里“重要性”理解为“便宜程度” # S1比S2同样便宜(1)S1比S3明显便宜(5)S2比S3明显便宜(5) matrix_c2 np.array([ [1, 1, 5], [1, 1, 5], [1/5, 1/5, 1] ]) # 对于准则C3交货期交货期越短越好 # S1比S2稍微快(3)S1比S3同样快(1)S2比S3稍微慢(1/3) matrix_c3 np.array([ [1, 3, 1], [1/3, 1, 1/3], [1, 3, 1] ]) criteria_matrices [matrix_c1, matrix_c2, matrix_c3] alternative_names [供应商S1, 供应商S2, 供应商S3] alternative_weights_matrix [] # 用于存储每个准则下各方案的权重 print(\n--- 方案层判断与权重计算 ---) for idx, (mat, criterion) in enumerate(zip(criteria_matrices, [质量, 价格, 交货期])): print(f\n对于准则 {criterion} 的判断矩阵:) print(mat) weights, lambda_max, ci, consistent, cr calculate_weights_and_consistency(mat) alternative_weights_matrix.append(weights) print(f 各方案权重: {weights}) print(f 一致性比率 CR: {cr:.4f}, 通过: {consistent}) # 3. 层次总排序计算每个方案的综合得分 alternative_weights_matrix np.array(alternative_weights_matrix).T # 转置使行对应方案列对应准则 total_scores np.dot(alternative_weights_matrix, criteria_weights) print(\n 层次总排序与决策结果 ) for name, score in zip(alternative_names, total_scores): print(f{name}: 综合得分 {score:.4f}) best_idx np.argmax(total_scores) print(f\n推荐选择{alternative_names[best_idx]} (最高综合得分: {total_scores[best_idx]:.4f})) if __name__ __main__: ahp_supplier_selection()运行这段代码你将看到完整的计算过程和结果输出。这个例子清晰地展示了AHP从准则权重计算到方案层权重合成最终得到总排序的整个流程。代码中包含了矩阵的构建、权重的计算、一致性检验以及结果展示形成了一个完整的闭环。实操心得在实际建模竞赛或应用中判断矩阵的获取往往不是通过程序交互输入而是来自问卷调研如专家打分。因此更常见的做法是将数据存储在Excel或CSV文件中然后使用pandas读取。我们可以很容易地修改get_judgment_matrix函数使其支持从DataFrame或文件路径加载矩阵这将大大提高代码的实用性。另外对于一致性检验不通过的情况除了提示警告更高级的实现可以尝试给出调整建议例如指出不一致性最可能源自哪几个判断值但这涉及到更复杂的算法如自动修正或灵敏度分析可以作为后续的扩展方向。4. TOPSIS法代码精讲从数据预处理到贴近度计算的标准化实现TOPSISTechnique for Order Preference by Similarity to Ideal Solution即“逼近理想解排序法”它的核心思想非常直观找到最优方案和最劣方案即理想解和负理想解然后通过计算每个方案与这两个参考点的距离来评估其优劣。与AHP的主观赋权不同TOPSIS常与客观赋权法如熵权法结合使用。代码实现的关键在于对数据预处理、距离计算和指标权重处理的精确把握。4.1 数据预处理归一化与指标类型的统一处理TOPSIS的第一步是将原始决策矩阵进行规范化以消除不同指标量纲和数量级的影响。最常用的方法是向量归一化Vector Normalization即每个元素除以该指标所有值的平方和的开方。此外我们必须明确每个指标是“效益型”越大越好还是“成本型”越小越好以便后续构造理想解。import numpy as np import pandas as pd from typing import List, Union def normalize_matrix(matrix: np.ndarray) - np.ndarray: 使用向量归一化方法处理决策矩阵。 参数: matrix: 原始决策矩阵形状为 (m个方案, n个指标)。 返回: 归一化后的决策矩阵。 # 向量归一化公式: zij xij / sqrt( sum(xij^2) ) # 对每一列每个指标进行计算 norm np.sqrt(np.sum(matrix ** 2, axis0)) # 避免除零错误如果某指标所有值都为0则归一化后全为0 norm[norm 0] 1 normalized_matrix matrix / norm return normalized_matrix def weighted_normalize_matrix(normalized_matrix: np.ndarray, weights: np.ndarray) - np.ndarray: 对归一化后的矩阵进行加权。 参数: normalized_matrix: 归一化后的决策矩阵。 weights: 权重向量长度等于指标数n。权重之和应为1。 返回: 加权规范化决策矩阵。 if not np.isclose(np.sum(weights), 1.0): print(f警告权重之和为{np.sum(weights)}将自动归一化。) weights weights / np.sum(weights) # 将权重向量扩展为与矩阵形状匹配然后进行元素乘法 # 这里使用 np.newaxis 来将一维权重向量转换为行向量便于广播 weighted_matrix normalized_matrix * weights[np.newaxis, :] return weighted_matrix这里有一个关键细节normalize_matrix函数中我们处理了norm 0的情况。这在什么情况下会发生当某个指标在所有方案上的值都完全相同时例如所有供应商的某个评分都是5其平方和开方后不为零。但当所有值都为0时就会出现除零错误。虽然在实际数据中不常见但健壮的代码必须考虑这种边界情况。接下来我们需要根据指标类型效益型/成本型来确定理想解和负理想解。我们设计一个函数来统一处理def identify_ideal_solutions(weighted_matrix: np.ndarray, indicator_types: List[str]) - Tuple[np.ndarray, np.ndarray]: 确定加权规范化矩阵的理想解和负理想解。 参数: weighted_matrix: 加权规范化决策矩阵形状 (m, n)。 indicator_types: 长度为n的列表每个元素为 benefit效益型或 cost成本型。 返回: ideal_best: 理想最优解正理想解形状 (n,)。 ideal_worst: 理想最劣解负理想解形状 (n,)。 n_indicators weighted_matrix.shape[1] if len(indicator_types) ! n_indicators: raise ValueError(指标类型列表长度必须与指标数量一致。) ideal_best np.zeros(n_indicators) ideal_worst np.zeros(n_indicators) for j in range(n_indicators): column weighted_matrix[:, j] if indicator_types[j] benefit: ideal_best[j] np.max(column) ideal_worst[j] np.min(column) elif indicator_types[j] cost: ideal_best[j] np.min(column) # 成本型指标最小值最优 ideal_worst[j] np.max(column) # 成本型指标最大值最劣 else: raise ValueError(f第{j}个指标类型必须为 benefit 或 cost 当前是 {indicator_types[j]}) return ideal_best, ideal_worst这个函数清晰地体现了TOPSIS的核心对于效益型指标理想解取最大值负理想解取最小值对于成本型指标则相反。明确区分指标类型是正确应用TOPSIS的前提很多初学者错误都源于此。4.2 距离计算与贴近度排序欧氏距离的几何意义计算每个方案到理想解和负理想解的距离通常采用欧几里得距离2-范数。然后计算每个方案与理想解的相对贴近度。def calculate_distances_and_closeness(weighted_matrix: np.ndarray, ideal_best: np.ndarray, ideal_worst: np.ndarray) - Tuple[np.ndarray, np.ndarray, np.ndarray]: 计算每个方案到理想解的距离、到负理想解的距离以及相对贴近度。 参数: weighted_matrix: 加权规范化决策矩阵。 ideal_best: 正理想解向量。 ideal_worst: 负理想解向量。 返回: dist_to_best: 每个方案到正理想解的距离形状 (m,)。 dist_to_worst: 每个方案到负理想解的距离形状 (m,)。 closeness: 每个方案的相对贴近度形状 (m,)。 m weighted_matrix.shape[0] dist_to_best np.zeros(m) dist_to_worst np.zeros(m) for i in range(m): # 计算方案i到正理想解的欧氏距离 dist_to_best[i] np.sqrt(np.sum((weighted_matrix[i, :] - ideal_best) ** 2)) # 计算方案i到负理想解的欧氏距离 dist_to_worst[i] np.sqrt(np.sum((weighted_matrix[i, :] - ideal_worst) ** 2)) # 计算相对贴近度 Ci dist_to_worst / (dist_to_best dist_to_worst) # 为了防止分母为零增加一个极小值 closeness dist_to_worst / (dist_to_best dist_to_worst 1e-12) return dist_to_best, dist_to_worst, closeness贴近度Ci的取值范围在0到1之间。Ci越大说明该方案离理想解越近离负理想解越远方案越优。这里我们给分母加了一个极小的数1e-12这是数值计算中防止除零错误的常用技巧。4.3 熵权法赋权让数据自己说话TOPSIS本身不产生权重权重需要外生给定。熵权法是一种客观赋权法它根据各指标数据所提供的信息量大小来确定权重。信息熵越小指标的变异程度越大提供的信息量越多其权重也应越大。def calculate_entropy_weights(decision_matrix: np.ndarray) - np.ndarray: 基于熵权法计算各指标权重。 参数: decision_matrix: 原始决策矩阵形状 (m, n)。要求所有元素为非负。 返回: 权重向量形状 (n,)。 # 1. 数据平移避免出现0值因为log(0)无定义。 # 通常平移一个非常小的数这里平移矩阵最小非零值的十分之一或直接平移1 matrix_shifted decision_matrix - np.min(decision_matrix) 1e-6 # 2. 计算每个方案在第j个指标下的比重 p_ij p matrix_shifted / np.sum(matrix_shifted, axis0, keepdimsTrue) # 3. 计算第j个指标的熵值 e_j # 当 p_ij 0 时根据极限p_ij * ln(p_ij) 0我们用np.where处理 with np.errstate(divideignore, invalidignore): entropy -np.sum(p * np.log(p), axis0) / np.log(decision_matrix.shape[0]) # 处理可能出现的NaN当某列所有p_ij都为0时 entropy np.nan_to_num(entropy, nan1.0) # 熵值最大为1表示完全无序无信息量 # 4. 计算差异系数 d_j 1 - e_j d 1 - entropy # 5. 计算权重 w_j d_j / sum(d_j) weights d / np.sum(d) return weights熵权法实现中有几个关键陷阱非负性要求熵权法要求决策矩阵元素非负。如果数据中有负数需要进行适当的平移处理如减去最小值使其全部变为非负。对数零问题计算熵值时需要取对数ln(p_ij)当p_ij0时无定义。虽然数学上规定0*ln(0) 0但代码中直接计算会得到NaN。我们的处理方法是先对矩阵进行微小平移1e-6确保所有p_ij 0。这是一种常用且稳定的做法。常数指标处理如果某个指标在所有方案上的值都完全相同变异度为0那么其熵值会达到最大值1或非常接近1差异系数d_j接近0从而导致其权重为0。这符合熵权法的逻辑不提供任何区分信息的指标权重应为0。代码中np.nan_to_num(entropy, nan1.0)将NaN替换为1正是为了处理这种极端情况。4.4 综合实战TOPSIS-熵权法评价城市发展水平让我们用一个具体的例子将数据预处理、熵权法赋权和TOPSIS排序串联起来。假设我们有4个城市方案从3个指标GDP增长率-效益型、失业率-成本型、PM2.5年均浓度-成本型来评价其发展水平。def topsis_city_evaluation(): 使用熵权法-TOPSIS综合评价城市发展水平。 print( 熵权法-TOPSIS城市发展水平评价 ) # 原始决策矩阵行-城市列-指标 # 指标顺序: [GDP增长率(%), 失业率(%), PM2.5浓度(μg/m³)] raw_data np.array([ [7.5, 3.8, 35], # 城市A [6.2, 5.1, 42], # 城市B [8.1, 4.2, 28], # 城市C [5.8, 6.0, 50], # 城市D ]) city_names [城市A, 城市B, 城市C, 城市D] indicator_names [GDP增长率, 失业率, PM2.5浓度] indicator_types [benefit, cost, cost] # 效益型成本型成本型 print(原始决策矩阵:) df_raw pd.DataFrame(raw_data, indexcity_names, columnsindicator_names) print(df_raw) # 步骤1: 熵权法计算权重 print(\n--- 步骤1: 熵权法计算指标权重 ---) entropy_weights calculate_entropy_weights(raw_data) for name, weight in zip(indicator_names, entropy_weights): print(f{name}: 权重 {weight:.4f}) # 步骤2: 数据归一化 normalized_matrix normalize_matrix(raw_data) print(\n--- 步骤2: 向量归一化后的矩阵 ---) df_norm pd.DataFrame(normalized_matrix, indexcity_names, columnsindicator_names) print(df_norm.round(4)) # 步骤3: 加权规范化 weighted_matrix weighted_normalize_matrix(normalized_matrix, entropy_weights) print(\n--- 步骤3: 加权规范化矩阵 ---) df_weighted pd.DataFrame(weighted_matrix, indexcity_names, columnsindicator_names) print(df_weighted.round(4)) # 步骤4: 确定理想解与负理想解 ideal_best, ideal_worst identify_ideal_solutions(weighted_matrix, indicator_types) print(f\n--- 步骤4: 理想解 (Z) ---) print(dict(zip(indicator_names, ideal_best.round(4)))) print(f--- 负理想解 (Z-) ---) print(dict(zip(indicator_names, ideal_worst.round(4)))) # 步骤5: 计算距离与贴近度 dist_best, dist_worst, closeness calculate_distances_and_closeness(weighted_matrix, ideal_best, ideal_worst) print(\n--- 步骤5: 距离计算与排序 ---) results [] for i, city in enumerate(city_names): results.append({ 城市: city, D (距理想解): dist_best[i], D- (距负理想解): dist_worst[i], 贴近度 Ci: closeness[i] }) results_df pd.DataFrame(results).set_index(城市) # 按贴近度降序排序 results_df results_df.sort_values(by贴近度 Ci, ascendingFalse) print(results_df.round(4)) print(f\n 综合评价排序结果 ) for rank, (city, row) in enumerate(results_df.iterrows(), start1): print(f第{rank}名: {city}, 贴近度 {row[贴近度 Ci]:.4f}) if __name__ __main__: topsis_city_evaluation()运行这段代码你将得到一个完整的TOPSIS分析报告。从输出中你可以清晰地看到熵权法根据数据本身的变异程度为三个指标赋予了不同的权重。通常数据离散程度越大的指标权重越高。加权规范化矩阵结合了指标权重是后续距离计算的基础。理想解和负理想解根据指标类型正确确定。最终根据贴近度Ci对城市进行排序Ci值最大的城市综合发展水平最优。注意事项与扩展在实际应用中原始数据可能包含极值或异常值这会对熵权法产生较大影响因为熵权法对数据分布敏感。因此在应用熵权法前进行适当的数据清洗和异常值处理是必要的。此外TOPSIS的距离计算除了欧氏距离也可以考虑使用曼哈顿距离等其他距离公式但欧氏距离是最为标准和常用的。如果指标间存在较强相关性可以考虑先使用主成分分析PCA等方法进行降维和去相关再用TOPSIS评价但这属于更高级的集成模型范畴了。这份代码提供了一个坚实、可复用的基础框架你可以根据具体问题的特点进行修改和扩展。