将一个数学优化模型翻译成Python语言，对其进行优化，并通过可视化解决方案来快速获取建模错误的反馈
👁️ 这是系列文章中的第3篇，涵盖了项目“Python中的智能决策支持系统应用于旅游”。我鼓励你看看它，以便对整个项目有一个整体的概述。如果你只对如何在Python中实现TSP模型感兴趣，那么你来对地方了：本文是独立的，我将带你逐步了解所有步骤——依赖安装、分析、代码、结果解释和模型调试。
本文延续了我们在上一篇文章中的旅程。在这里，我们采用前一篇文章中制定的数学模型，并使用Pyomo在Python中实现它，遵循良好的实践。然后，优化该模型，并可视化和分析其解决方案。出于教育目的，我们发现初始模型制定不完整，所以我将展示如何从基本原理推导出修复制定所需的约束条件。这些新的约束条件添加到Pyomo模型中，再次进行分析和验证新的解。
目录
1. 使用Pyomo在Python中实现模型
- 1.1. 安装依赖项
- 1.2. 数学变成代码
2. 解决和验证模型
- 2.1. 解决模型
- 2.2. 可视化结果
- 2.3. 分析结果
3. 修复制定
- 3.1. 动机思路
- 3.2. 将动机思路表达为逻辑蕴含
- 3.3. 将逻辑蕴含形式化为线性约束
- 3.4. 推导“大M”的适当值
4. 实施和验证新制定
- 4.1. 增强Pyomo模型
- 4.2. 绘制更新后模型的解决方案
5. 结论（下一个迭代中）
📖 如果你没有读过以前的文章，别担心。数学制定也在这里声明（但没有推导），每个模型组件都与其代码实现相邻。如果你不理解事物的来源或含义，请阅读“迭代2”的文章，如果需要更多上下文，请阅读“迭代1”的文章。
1. 使用Pyomo在Python中实现模型
Python被用作数据科学中的顶级语言，而Pyomo则是处理大规模模型效果最好的开源库之一。
在这一部分中，我将逐个介绍公式中定义的模型组成部分，并解释如何将其转换为Pyomo代码。我尽量不留下任何空白，但如果你有其他意见，请在评论中提出问题。
免责声明：目标读者预计对Pyomo甚至建模都是新手，因此重视简明和简单的实现，而非编程最佳实践。现在的目标是教授优化建模，而非软件工程。随着这个概念验证进一步发展成为一个更复杂的MVP，代码将在未来迭代中逐步改进。
1.1 安装依赖项
对于着急的人
安装（或确保已经安装了）库pyomo、networkx和pandas以及包glpk。
📝 包含GLPK求解器的
glpk包是一个（开源的）外部求解器，我们将用它来优化我们创建的模型。Pyomo用于创建问题模型并将其传递给GLPK，后者将运行算法进行优化处理。然后，GLPK将解返回给Pyomo模型对象，并将其存储为模型属性，因此我们可以方便地在Python中使用它们而不离开。
推荐的安装GLPK的方法是通过conda，这样Pyomo可以轻松找到GLPK求解器。要一次性安装所有依赖项，请运行：
conda install -y -c conda-forge pyomo=6.5.0 pandas networkx glpk对于有条理的人
我建议创建一个单独的虚拟环境，在其中安装所有需要遵循本系列文章的库。复制以下文本：
name: ttp # 旅行旅游问题channels: - conda-forgedependencies: - python=3.9 - pyomo=6.5.0 - pandas - networkx - glpk # 用于优化模型的外部求解器 - jupyterlab # 如果你不使用Jupyter Lab作为IDE，则注释掉这一行并将其保存为名为environment.yml的YAML文件。在相同位置打开Anaconda命令提示符，然后运行以下命令：
conda env create --file environment.yml几分钟后，环境将被创建，并安装所有依赖项。运行conda activate ttp进入该环境，在终端中启动Jupyter Lab（通过运行jupyter lab），开始编码！
1.2. 数学变成代码
首先，确保Pyomo可以找到GLPK求解器
### ===== 代码块 3.1 ===== ###import pandas as pdimport pyomo.environ as pyoimport pyomo.versionimport networkx as nxsolver = pyo.SolverFactory("glpk")solver_available = solver.available(exception_flag=False)print(f"求解器'{solver.name}'可用：{solver_available}")if solver_available: print(f"求解器版本：{solver.version()}")print("pyomo版本：", pyomo.version.__version__)print("networkx版本：", nx.__version__)
求解器'glpk'可用：True求解器版本：(5, 0, 0, 0)pyomo版本：6.5networkx版本：2.8.4⛔ 如果你收到了消息
'glpk' available: False，表示求解器没有正确安装。请尝试以下措施解决问题：– 仔细重新执行安装步骤
– 在基本（默认）环境中运行
conda install -y -c conda-forge glpk– 尝试安装适合你的其他求解器
然后读取距离数据的CSV文件
### ===== 代码块 3.2 ===== ###path_data = ( "https://raw.githubusercontent.com/carlosjuribe/" "traveling-tourist-problem/main/" "Paris_sites_spherical_distance_matrix.csv")df_distances = pd.read_csv(path_data, index_col="site")df_distances
现在我们进入[敏捷运筹学工作流程]的“第四阶段”，如下图所示：
任务是根据之前创建的数学模型，将其按照数学定义精确地实现为代码。
👁️ 如果这样做会使模型实现更容易，我们可以创建任意数量的Python对象，但在编码过程中不能修改基本模型。 这会导致数学模型和计算机模型不同步，从而使后续的模型调试变得非常困难。
我们实例化一个空的Pyomo模型，将模型组件存储为属性：
model = pyo.ConcreteModel("TSP")1.2.1. 集合
为了创建站点集𝕊 = {卢浮宫，埃菲尔铁塔，…，酒店}，我们从数据帧的索引中提取它们的名称，并用它创建一个名为sites的Pyomo Set：
### ===== 代码块 3.3 ===== ###list_of_sites = df_distances.index.tolist()model.sites = pyo.Set(initialize=list_of_sites, domain=pyo.Any, doc="要访问的所有站点的集合（𝕊）")为了创建派生集合
我们将筛选器 𝑖 ≠ 𝑗 存储在一个构造规则中（Python函数_rule_domain_arcs），并将此规则传递给初始化Set时的filter关键字。请注意，此筛选器将应用于站点的笛卡尔积（𝕊 × 𝕊），并过滤掉不满足规则的笛卡尔积成员。
### ===== 代码块 3.4 ===== ###def _rule_domain_arcs(model, i, j): """连接站点（𝔸）的所有可能弧""" # 仅当站点i和站点j不同时才创建对(i, j) return (i, j) if i != j else None rule = _rule_domain_arcsmodel.valid_arcs = pyo.Set( initialize=model.sites * model.sites, # 𝕊 × 𝕊 filter=rule, doc=rule.__doc__)1.2.2. 参数
参数
𝐷ᵢⱼ ≔ 站点𝑖和站点𝑗之间的距离
使用构造函数pyo.Param创建，该构造函数以第一个（位置）参数为索引它的域𝔸（model.valid_arcs），关键字参数initialize为另一个构造规则_rule_distance_between_sites，对于𝔸中的每对（𝑖，𝑗），都会对其进行评估。 在每次评估中，从数据帧df_distances中获取距离的数值，并将其与对（𝑖，𝑗）进行内部关联：
### ===== Code block 3.5 ===== ### def _rule_distance_between_sites(model, i, j): """ site i和site j之间的距离（𝐷𝑖𝑗）""" return df_distances.at[i, j] # 从数据框中获取距离 rule = _rule_distance_between_sitesmodel.distance_ij = pyo.Param(model.valid_arcs, initialize=rule, doc=rule.__doc__)1.2.3. 决策变量
由于ηᵢⱼ具有与𝐷ᵢⱼ相同的“索引域”，构建该组件的方式非常相似，唯一的区别是在此处不需要构造规则。
model.delta_ij = pyo.Var(model.valid_arcs, within=pyo.Binary, doc="是否从site i到site j (𝛿𝑖𝑗)")pyo.Var的第一个位置参数保留给索引集𝔸, 并且变量的"类型"由关键字参数within指定。在这里，𝛿ᵢⱼ的取值范围只是0和1，所以它是二元的。数学上，我们会写成𝛿ᵢⱼ ∈ {0, 1}，但是我们可以直接在Pyomo中通过在创建变量时设置
within=pyo.Binary来表达这一点，而不是创建单独的约束来指示这一点。1.2.4. 目标函数
要构建目标函数，我们可以用一个函数“存储”表达式，该函数将用作目标的构建规则。这个函数只有一个参数，模型，用于获取构建表达式所需的任何模型组件。
### ===== Code block 3.6 ===== ### def _rule_total_distance_traveled(model): """ 总旅行距离 """ return pyo.summation(model.distance_ij, model.delta_ij) rule = _rule_total_distance_traveled model.obj_total_distance = pyo.Objective(rule=rule, sense=pyo.minimize, doc=rule.__doc__)观察数学表达式和函数的返回语句之间的对应关系。我们使用
sense关键字指定这是一个最小化问题。1.2.5. 约束条件
如果您回忆上一篇文章，一种方便的方法是通过同时执行每个站点“进入”一次和“离开”一次，来确保每个站点仅访问一次。
每个站点只能进入一次
def _rule_site_is_entered_once(model, j): """ 每个站点j必须由另一个站点访问一次 """ return sum(model.delta_ij[i, j] for i in model.sites if i != j) == 1 rule = _rule_site_is_entered_oncemodel.constr_each_site_is_entered_once = pyo.Constraint( model.sites, rule=rule, doc=rule.__doc__)每个站点只能离开一次
def _rule_site_is_exited_once(model, i): """ 每个站点 i 必须离开到正好一个其它站点 """ return sum(model.delta_ij[i, j] for j in model.sites if j != i) == 1rule = _rule_site_is_exited_oncemodel.constr_each_site_is_exited_once = pyo.Constraint( model.sites, rule=rule, doc=rule.__doc__)1.2.6. 模型的最终检查
模型的实现已经完成。要查看整个模型的样子，我们应该执行
model.pprint()，并且四处浏览一下，看是否遗漏了一些声明或者引起了一些错误。为了了解模型的规模，我们打印出它包含的变量和约束的数量：
def print_model_info(model): print(f"名称：{model.name}", f"变量数量：{model.nvariables()}", f"约束数量：{model.nconstraints()}", sep="\n- ")print_model_info(model)#[输出]:# 名称：TSP# - 变量数量：72# - 约束数量：18约束或变量少于100个，这是一个小规模的问题，求解器将会相对快速地进行优化。
2. 求解和验证模型
2.1. 求解模型
[AOR流程图]中的下一步是对模型进行优化，并检查解决方案：
res = solver.solve(model) # 优化模型print(f"找到最优解：{pyo.check_optimal_termination(res)}")# [输出]: 找到最优解：真好消息！求解器已经找到了这个问题的一个最优解！让我们来检查一下，这样我们在到达巴黎时就知道应该遵循哪条线路！
为了进行一个非常快速的检查，我们可以运行
model.delta_ij.pprint()，它将打印出所有（最优的）𝛿ᵢⱼ变量的值，值要么是0，要么是1：仅仅通过查看已选择的弧线列表来可视化一条路径是困难的，更不用说分析它以验证我们是否正确地建模了。
要真正理解解决方案，我们需要将其可视化。
2.2. 可视化结果
一图胜千言
由于我们涉及到节点和弧线，最简单的可视化解决方案的方式是将其绘制为图形。请记住，这是一个概念验证，因此快速有效的反馈胜过美观。更深刻的可视化可以等到我们有一个完善的MVP。现在，让我们编写一些辅助函数来高效绘制解决方案。这
函数
extract_solution_as_arcs接受一个已解决的Pyomo模型，并从解决方案中提取出”已选择的弧线”。接下来，函数plot_arcs_as_graph将活跃弧线的列表存储在一个图形对象中，以便更容易进行分析，并绘制该图形，以使酒店成为唯一的红色节点，作为参考。最后，函数plot_solution_as_graph调用上述两个函数，将给定模型的解决方案显示为图形。### ===== 代码块 3.7 ===== ### def extract_solution_as_arcs(model): """ 从解决的模型中提取活动（选定的）弧的列表， 只保留其二进制变量为1的弧""" active_arcs = [(i, j) for i, j in model.valid_arcs if model.delta_ij[i, j].value == 1] return active_arcs def plot_arcs_as_graph(tour_as_arcs): """ 接受一个表示弧的元组列表，将其转换为网络图并绘制""" G = nx.DiGraph() G.add_edges_from(tour_as_arcs) # 将解决方案存储为图形 node_colors = ['red' if node == 'hotel' else 'skyblue' for node in G.nodes()] nx.draw(G, node_color=node_colors, with_labels=True, font_size=6, node_shape='o', arrowsize=5, style='solid') def plot_solution_as_graph(model): """ 将给定模型的解决方案绘制为图形""" print(f"总距离：{model.obj_total_distance()}") active_arcs = extract_solution_as_arcs(model) plot_arcs_as_graph(active_arcs)现在我们可以看到解决方案的真正样子：
plot_solution_as_graph(model)显然，这不是我们的预期结果！没有一个单一的旅行途经所有景点并返回酒店。确实，所有景点都被访问了，但只作为小的离散的景点群的一部分。从技术上讲，我们指定的约束性被仔细遵守：每个景点只进入一次和退出一次，但是整体结果并不是我们预期的一个单一的旅行，而是一组子旅行。这意味着我们在之前的文章中做出的假设是错误的，还有一些其他模型缺少的东西，将“子旅行不允许在解决方案中”编码。
2.3. 分析结果
出了什么问题？
当模型的解决方案没有意义时，仅有一种可能的解释：模型是错误的¹。
求解器给我们了一个真正的最优解决方案，但我们给它一个与我们要解决的问题不匹配的模型。现在的任务是找出为什么，以及我们犯了什么错误的地方。经过反思，显而易见的候选人是我们在上一篇文章的第4.4节的最后两段中所做的可疑假设，我们在其中设计了数学模型。我们（现在我们知道是错的）假设一个单一旅行的形成将自然而然地遵循确保每个景点仅访问一次的两个约束条件。但正如我们刚刚可视化的那样，它没有。为什么？
错误的根本原因在于我所称的“不言而喻的常识”：我们对世界的了解是如此明显，以至于我们忘记在模型中明确说明它
我们知道，暗示地，在参观景点时是不可能进行瞬时移动的，但我们没有明确告诉模型。这就是为什么我们观察到这些小的子旅行，连接其中的一些景点，但并非全部。模型“认为”从一个景点瞬时传送到另一个景点是可以的，只要它在一个景点停留一次并进入一次（再次参见图3.2）。如果我们看到子旅行，那只是因为我们告诉模型要最小化旅行的距离，而瞬时传送在节省距离方面很有帮助。
因此，我们需要防止这些子路径的形成，以获得一个现实的解决方案。我们需要设计一些新的约束条件，“告诉模型”子路径是禁止的，或者等效地说，解决方案必须是一条单一的路径。让我们选择后者，并从第一原理推导出一组直观且有效的约束条件。
3. 修正公式
参考[敏捷运营研究流程]，我们现在处于模型重新制定阶段。重新制定模型可能是为了改进或修复模型。我们的目标是将解决方案强制转化为一条从起始点（酒店）开始并结束于酒店的单一路径。挑战在于如何将该要求编码为一组线性约束条件。下面是一个基于路径特性的想法。
3.1. 动因理念
我们有𝑁个要“遍历”的地点，包括酒店。因为我们从酒店开始，这意味着还剩下𝑁-1个要访问的地点。如果我们追踪我们访问这些地点的“时间顺序”，即第一个目的地（酒店之后）编号为1，第二个目的地编号为2，依此类推，那么在返回酒店之前的最后一个目的地将被编号为𝑁-1。如果我们称这些用于追踪访问顺序的编号为“等级”，那么在路径中，任何地点（除酒店之外）的等级始终比其前面的地点的等级高1个单位。如果我们能制定一组约束条件，使得任何可行解都满足这种模式，我们将引入一种“时间顺序”要求。事实证明我们确实可以做到。
3.2. 将动因理念表达为逻辑蕴含
💡 这就是我们要求任何可行解满足的“模式”：
任何地点（酒店以外的地点）的等级必须始终比其前面的地点的等级高1个单位。
我们可以将这个模式重新表达为逻辑蕴含，如下：“如果且仅当𝑗紧邻𝑖被访问时，地点𝑗的等级必须比地点𝑖的等级高1个单位，对于所有不包括酒店𝐻的弧(𝑖,𝑗)”。这个措辞在数学上被表达为：
其中𝑟ᵢ是我们需要跟踪的（尚未知的）访问顺序的新变量。为了与决策变量区分开来，我们将其称为“等级”变量。右侧的读法是“对于所有𝑖和𝑗属于除酒店之外的所有地点集合”。为了符号的方便，我们定义新集合𝕊*来存储除酒店（用𝐻表示）之外的所有地点：
这使我们能够简洁地定义等级变量为：
👁️酒店没有相关的等级变量，因为它将同时成为任何旅行的起点和终点，这违反了“旅行中等级变量递增”的模式。因此，每个网站的等级始终被强制在采取新的弧时增加，从而确保除非环路在唯一没有等级变量的网站（即酒店）关闭，否则不允许封闭的循环
𝑟ᵢ的边界是从其描述中得出的：等级从1开始，并且单调增加，直到访问𝕊*中的所有网站，最后停止在| 𝕊* |（非酒店网站集合的大小）上。此外，我们允许它们取任何正实数值：
3.3. 将逻辑蕴含式转化为线性约束条件
现在的挑战是将这个逻辑蕴含式转化为一组线性约束条件。谢天谢地，线性不等式也可以用于强制逻辑蕴含关系，而不仅仅是有限资源的限制。
一种方法是通过所谓的大M法，通过这种方法声明约束条件，当满足您关心的条件时，约束条件适用（变为活动状态），而不满足您关心的条件时，约束条件变得多余（变为非活动状态）。该技术之所以被称为“大M”，是因为它使用常数𝑀，而𝑀足够大，以至于当出现在约束条件中时，在每种情况下都使其多余。当约束条件中不出现𝑀时，约束条件是“活动的”，即它在强制所需的蕴含关系。
但是，是什么决定了约束条件是否“活动”？简短的答案是约束条件适用于决策变量本身的值。让我们看看它是如何工作的。
所需的蕴含关系是当 𝛿ᵢⱼ = 1 时，𝑟ⱼ = 𝑟ᵢ + 1。我们可以将表达式中的 1 替换为 𝛿ᵢⱼ，得到
这是我们希望当 𝛿ᵢⱼ = 1 时成立的关系，但是当 𝛿ᵢⱼ = 0 时不成立。为了“纠正”当 𝛿ᵢⱼ = 0 时的情况，我们添加了一个冗余项，𝑅𝑇，其作用是仅在 𝛿ᵢⱼ = 0 时“停用约束条件”。因此，这个冗余项必须包含变量𝛿ᵢⱼ，因为它依赖于它。
在这个背景下，“停用约束条件”意味着“使其冗余”，因为冗余约束条件在模型中与不存在的约束条件具有相同的效果。
让我们来看一下如何推导出RT的表达式。RT的表达式需要满足以下属性：
满足条件（1）需要𝑅𝑇(𝛿ᵢⱼ = 1) = 0，因此𝑅𝑇的表达式必须包含乘法因子(1-𝛿ᵢⱼ)，因为当𝛿ᵢⱼ = 1时，它变为0。这个形式使得当𝛿ᵢⱼ = 1时，𝑅𝑇“消失”，或者当𝛿ᵢⱼ = 0时，“被减少”为一个常数（我们称之为M）。因此，冗余术语的一个候选是
其中𝑀应根据问题数据确定（稍后详细说明）。
为了满足条件（2）对于所有可能的𝑖和𝑗，我们需要将表达式（3.11）中的等式变为不等式（=变为≥），并找到一个充足大的常数𝑀的值，以便无论𝑟ⱼ和𝑟ᵢ取什么值，约束条件始终得到满足。这就是“大M”中的“大”的含义。
一旦我们找到这样一个足够大的常数𝑀，我们的“逻辑蕴含”约束条件将采取以下形式
将这些约束条件引入模型将会强制求解结果为单一路径。但是，如果我们不先指定M的一个好值，约束条件将不会产生期望的效果。
3.4. 推导“大M”的适当值
由于目标是使𝑅𝑇(𝛿ᵢⱼ = 0) = 𝑀，我们可以通过将𝛿ᵢⱼ = 0代入表达式（3.12）中的一般约束条件来推导出𝑀的适当值：
为了对于所有非酒店网站𝑖，𝑗，𝑟ⱼ − 𝑟ᵢ ≥ 𝑀能够满足，我们需要𝑟ⱼ − 𝑟ᵢ的下界也大于𝑀。𝑟ⱼ − 𝑟ᵢ的下界（LB）是𝑟ⱼ − 𝑟ᵢ可能取的最小值，可以通过以下方式获得：
其中 𝑟ᵐⁱⁿ 是最小可能等级，而 𝑟ᵐᵃˣ 是最大可能等级。因此，为了使不等式（1）在所有站点的所有等级成立，以下不等式也必须成立：
借助这个不等式，我们知道 𝑀 必须取的最小值，以使大-M方法有效，即
那么 𝑟ᵐⁱⁿ 和 𝑟ᵐᵃˣ 的值是多少呢？根据我们的约定，我们将第一个访问的站点的等级设为1，这当然是最小等级（即 𝑟ᵐⁱⁿ = 1）。由于每访问一个站点等级增加1个单位，游览中的最后一个非酒店站点将具有最大等级，等于所有非酒店站点的总数。由于非酒店站点的数量可能不同，我们需要一个一般的表达式。如果你还记得的话，我们定义了集合 𝕊* 包含所有非酒店站点，所以我们需要的数量是集合 𝕊* 的大小（即元素的数量），在数学符号中用 | 𝕊* | 表示。因此，我们推断出 𝑟ᵐⁱⁿ = 1 且 𝑟ᵐᵃˣ = | 𝕊* |。将其代入到Expression (3.14)中，我们最终得到了𝑀的适当值：
由于 𝕊* 总是有超过2个待访问站点（否则，一开始就没有决策问题），所以在这个模型中，“big M”值总是一个负的“大”数值。
📝 理论值与计算值对比
理论上，我们可以任意选择比这里推导出的 𝑀 更“负”的值，甚至编造巨大的数以确保它们足够大，以避免此计算，但这不是好的实践。如果 𝑀 变得太大（绝对值上），它可能会在求解器算法中产生性能问题，或者在最坏的情况下，甚至使求解器将不可行的解视为可行解。因此，推荐的做法是从问题的数据中推导出一个紧致但足够大的值作为 𝑀。
现在，我们已经推导出了“big M”的适当值，我们将其存储在一个新的模型参数中，以便更容易重用。有了这个，解决子旅游路径消除约束的集合已经准备好了，它的“完整形式”如下：
为了保持事物的透视，注意这实际上是之前我们在表达式(3.6)中制定的原始“约束等价物”的表示：
恭喜！我们终于有了一组可以添加到我们的模型中的约束条件。想出这些约束条件是困难的部分。现在，让我们验证将它们添加到模型中确实会导致旅行子路线的消失。
4. 实施并验证新的公式
4.1. 使用新公式增强Pyomo模型
修订和更正模型需要添加一些更多的集合、参数、变量和约束条件。让我们按照初始制定阶段的相同顺序将这些新的模型组件添加到Pyomo模型中。
4.1.1. 集合和参数
- 感兴趣的地点，𝕊*：所有地点的集合不包括酒店：
Pyomo Set对象具有与Python集合兼容的操作，因此我们可以直接在Pyomo集合和Python集合之间进行差异操作：
model.sites_except_hotel = pyo.Set( initialize=model.sites - {'hotel'}, domain=model.sites, doc="感兴趣的地点，即除了酒店之外的所有地点（𝕊*）")- 大M，用于旅行子路线消除约束的新参数：
model.M = pyo.Param(initialize=1 - len(model.sites_except_hotel), doc="用于使某些约束多余的大M")4.1.2. 辅助变量
- 排名变量，rᵢ：用于跟踪访问地点的顺序：
model.rank_i = pyo.Var( model.sites_except_hotel, # i ∈ 𝕊*（索引） within=pyo.NonNegativeReals, # rᵢ ∈ ℝ₊（域） bounds=(1, len(model.sites_except_hotel)), # 1 ≤ rᵢ ≤ |𝕊*| doc="每个地点的排名，用于跟踪访问顺序")在注释中，您可以看到完整数学定义的各个元素如何与Pyomo变量声明函数pyo.Var的参数对应得很好。我希望这可以帮助您在开始构建Pyomo模型之前理解拥有一个定义明确的数学模型的价值。实施将会更自然流畅，并且错误可能性更小。
4.1.3. 约束条件
- 解决方案必须是在酒店开始和结束的单一旅行：
def _rule_path_is_single_tour(model, i, j): """对于每一对非酒店站点(i, j)，如果站点j是从站点i访问的，则j的排名必须严格大于i的排名。""" if i == j: # 如果站点重合，则跳过创建约束 return pyo.Constraint.Skip r_i = model.rank_i[i] r_j = model.rank_i[j] delta_ij = model.delta_ij[i, j] return r_j >= r_i + delta_ij + (1 - delta_ij) * model.M# 非酒店站点的交叉乘积，用于索引约束non_hotel_site_pairs = model.sites_except_hotel * model.sites_except_hotelrule = _rule_path_is_single_tourmodel.constr_path_is_single_tour = pyo.Constraint( non_hotel_site_pairs, rule=rule, doc=rule.__doc__)Pyomo模型已更新。在添加子路径消除约束后，它增长了多少？
print_model_info(model)# [返回]:# 名称：TSP# - 变量数量：80# - 约束数量：74我们从72个变量增加到80个变量，从18个约束增加到74个约束。显然，这个公式中约束的数量比变量多，因为它将我们之前的约束数量增加了四倍。这通常是使模型更“现实”的“代价”，因为从实际角度看，通常需要限制可行解的数量（如果数据不变）。
我们可以像往常一样使用model.pprint()来检查模型结构。尽管随着模型组件数量的增长，这种“打印”很快就会变得毫无价值，但它仍然可以让我们快速了解模型的总体情况和规模。
4.2. 绘制更新后模型的解决方案
让我们解决更新后的模型并绘制新的解决方案。祈祷。
res = solver.solve(model) # 优化模型solution_found = pyo.check_optimal_termination(res)print(f"找到最优解：{solution_found}")# [返回]: 找到最优解：Trueif solution_found: plot_solution_as_graph(model)
现在让我们来谈谈！这正是我们在添加新约束之后所期望的：没有形成任何子路径，使得解决方案变成了一个单一的旅游线路。
请注意，显然，这个求解模型的目标值现在是14.9公里，而不是我们在不完整的模型中得到的虚假的5.7公里。
👁️ 图形的绘制并不是地图上的路径
请注意，此图像仅是图形的一个可能的绘制，而不是地理轨迹。您在图中看到的圆圈和链接不对应地理空间中的实际路径（如果我们没有在其创建中使用任何地理信息，那它又如何可能呢？）。您可以通过多次运行
plot_solution_as_graph(model)来验证这一点：每次运行它时，节点的位置都会发生变化。图是抽象的数学结构，通过“链接”将“点”连接在一起，表示任何种类的实体之间的任何种类的关系。我们在这里使用了一个图来研究解决方案的有效性，而不是在巴黎展示真实旅游线路。我们在[这篇文章]中进行。
5. 结论（或下一步的规划）
通过对解决方案进行最后验证，我们可以得出结论：这个更新的模型可以解决旅行商问题的任何实例，因此我们可以将其视为成功的概念验证（POC）。
💡 一步一步演变的解决方案
这个POC尚不能解决我们最初的复杂旅游问题，但它确实解决了我们作为第一个阶段向更复杂解决方案迈出的最小有价值问题。因此，它使我们从循证（基于证据）的角度更接近解决更复杂问题的有价值解决方案。有了一个最小的工作示例，我们可以更好地评估在我们想要的方向前进时需要做什么，以及哪些临时简化方法可以被暂时地舍弃，直到该解决方案的更成熟版本到来。在任何时候都有一些有用的东西，我们将开发一个越来越有价值的系统，直到满意为止。这就是“敏捷之路”有效解决方案的本质。
通过这种方法的有效性得到证明，我们必须扩展和改进它，逐渐包含更多原始问题的特性，每次迭代提供越来越有价值的解决方案。在这个 POC 中，我们专注于设计和构建一个基本模型，所以我们必须假设一组固定的站点和它们的距离矩阵作为给定。当然，这是有限制的，下一步应该是有一个模型，可以接受任意数量的站点。为此，我们需要一种自动根据站点列表和它们的地理坐标生成距离矩阵的方法。这是我们下一次冲刺的目标。
5.1. 接下来是什么
在下一篇文章（冲刺4）中，我们将开发一个类，可以自动从任何站点列表生成距离矩阵。当与我们刚刚构建的模型结合使用时，这个功能将允许我们在其他方面快速解决不同输入的多个模型并进行比较。此外，以这种方式推广解决方案将在未来的冲刺中使我们的工作更加轻松，当我们进行一些敏感性和场景分析时，这将非常有帮助。另外，我们将将这个概念验证升级为“MVP状态”，我们将开始使用面向对象的代码来保持良好的组织性，并为扩展性做好准备。不要错过机会，立即开始吧！
脚注
- 实际上，出错的原因还有另一个：模型所依赖的数据也可能存在错误，而不仅仅是模型的构建方式。但在某种程度上，如果将“模型”看作是一个“模型实例”，即具体的模型与具体的数据，那么如果数据有误，模型当然也会出错，这正是我此前的意图。
感谢阅读，期待下次再见！📈😊
请随时关注我，向我提问，在评论中给我反馈，或通过LinkedIn联系我。