当前位置：首页 > news >正文

数控切割路径优化：双层旅行商问题与迭代求解策略

news 2026/8/26 23:01:57

1. 问题引入：从一块钢板到最优切割路径

五一数学建模竞赛的A题，每年都像一道精心设计的“工业谜题”，今年这道“钢板最优切割路径问题”也不例外。乍一看，题目描述的是数控切割机如何在一块钢板上高效地切割出多个零件，目标是找到一条让空程（切割头不进行切割的移动距离）最短的路径。这听起来像是一个经典的“旅行商问题”（TSP）变种——切割头需要“访问”每一个待切割的图形轮廓，但又不完全一样，因为切割头在每个图形内部还需要沿着轮廓走完一圈。很多初次接触的同学可能会直接套用现成的智能优化算法，比如遗传算法或者模拟退火，去优化图形的访问顺序。但如果你真这么做了，很可能第一步就走偏了。这道题的核心难点和魅力，恰恰在于它比单纯的TSP多了一层：图形内部的切割起点选择。

想象一下，你手里拿着一支焊枪（切割头），面前是一块大钢板，上面用粉笔画好了十几个不同形状的零件轮廓。你的任务是把它们都切下来。你当然可以决定先切圆，再切方，最后切那个复杂的五角星（这就是访问顺序）。但当你移动到圆旁边准备开切时，一个新的问题出现了：从圆周上的哪一点开始下刀？你从圆的正上方点开始顺时针切，和从正下方点开始逆时针切，对于完成这个圆本身来说没有区别，但对于你切完这个圆后，前往下一个方形的“空程”距离，影响可就大了。这个起始点，我们称之为“切入点”或“切割起点”。最优路径 = 最优的外部图形访问序列 × 每个图形内部的最优切割起点。这是一个典型的双层优化问题，两层决策相互耦合，这才是题目真正的“骨头”。

我见过很多解题思路，要么只优化序列，默认从图形上某个固定点（比如离上一个图形最近的点）开始切；要么试图用超级复杂的算法同时优化两者，结果陷入组合爆炸，算力撑不住。今天，我想分享一套经过实战检验的、清晰且可实现的思路：“分而治之，迭代逼近”。我们不追求一步到位的全局最优，而是通过合理的建模和高效的局部搜索，找到一个质量非常高、逻辑清晰的可行解。下面，我就把这套方法的每一个环节掰开揉碎，包括其背后的原理、具体的实现步骤、我踩过的坑，以及能直接运行的Python参考代码。

2. 模型构建：如何将钢板切割抽象为数学模型

面对一个实际问题，首要任务是进行合理的抽象和简化，建立数学模型。我们不能一头扎进代码里，必须先想清楚用什么数据结构来表示问题，以及要优化的是什么。

2.1 关键概念定义与假设

首先，我们明确几个关键概念和必要的假设，这能让模型更清晰，也便于后续编程：

切割头：视为一个没有大小的点。这在路径规划中是常见假设，简化了碰撞检测（本题未强调碰撞，故更可忽略）。
空程：切割头在非切割状态下的移动距离。核心优化目标就是最小化总空程。切割路径本身的长度是固定的（即所有零件轮廓周长之和），无法优化，因此我们只关心空程。
零件图形：题目中通常是多边形（矩形、三角形等）和圆形。我们需要用数学方式描述它们。
- 多边形：用一组有序的顶点坐标表示，例如矩形[(x1,y1), (x2,y2), (x3,y3), (x4,y4)]。切割路径就是依次连接这些顶点并回到起点。
- 圆形：用圆心(cx, cy)和半径r表示。切割路径是一个圆。为了统一处理，我们需要将圆形离散化，即用正多边形来近似。例如，用360个点来近似一个圆，足以保证精度。
切入点：在每个图形轮廓上开始切割的点。对于多边形，可以是任意顶点；对于圆形，是离散化后多边形上的任意点。
切割方向：对于多边形，通常假设固定为顶点给定的顺序（如顺时针）。改变顺序可能意味着“翻转”零件，在实际切割中可能不被允许，因此我们默认固定。对于圆形，方向（顺时针/逆时针）不影响空程。

基于以上，我们可以做出一个至关重要的建模决策：将每个零件的完整闭合轮廓的切割，转化为对一个特殊的“城市”的访问。这个“城市”不是一个点，而是一条闭合的环。访问这个“城市”意味着：1. 移动切割头到该环上的某个点（切入点）；2. 完整遍历这个环（执行切割）；3. 从环上的终点（也就是切入点，因为闭合）离开，前往下一个“城市”。

2.2 数学模型的形式化表述

设共有N个零件。对于零件i：

其轮廓离散化或表示为一系列有序点P_i = [p_i1, p_i2, ..., p_iM_i]，其中p_i1和p_iM_i首尾相连。M_i是描述该轮廓的点数。
定义d_cut(i)为切割零件i所需的固定路径长度，即其轮廓周长。
定义enter_point(i, k)为选择零件i的第k个点 (p_ik) 作为切入点。

我们需要决策：

一个零件的排列顺序π = [π1, π2, ..., πN]，表示切割的先后次序。
对于每个零件πj，选择一个切入点索引k_j。

目标函数：最小化总空程。总路径 = 初始空程(从原点或初始点到第一个零件切入点) + Σ(零件间空程) + Σ(零件切割固定长度)。由于 Σ(零件切割固定长度) 是常数，因此优化目标等价于：Minimize: D_start(enter(π1)) + Σ D(exit(πj), enter(π(j+1)))其中：

D_start(p)是从切割头初始位置到切入点p的距离。
exit(πj)是零件πj切割结束的点，由于切割是闭合的，它等于该零件的切入点enter(πj)。
D(a, b)是点a到点b的欧几里得距离。

至此，我们成功地将一个复杂的物理切割问题，转化为了一个清晰的组合优化问题：为一个特殊的“旅行商”规划路线，每个“城市”允许你选择其边界上的任何一个点作为“访问站”。

3. 核心求解策略：分治与迭代逼近框架

直接求解上述双层优化问题非常困难。我们的策略是将其分解，并迭代改进。

3.1 策略一：固定切入点，优化访问序列

这是简化问题的第一步。我们暂时“冻结”内部优化，为每个零件i预先选定一个切入点。一个直观简单的策略是：选择每个零件轮廓上距离钢板中心最近的点，或者距离上一个零件“可能位置”最近的点。但在迭代初期，我们没有序列信息，一个稳健的初始选择是：选择每个零件轮廓上，距离所有其他零件轮廓最近点的平均位置最近的那个点。简单起见，初始阶段我们可以选择每个零件的一个特征点（如多边形重心、圆心）作为其“代表点”，但注意，空程计算时，移动是从一个零件的切入点到另一个零件的切入点，而不是重心到重心。

假设我们已经为每个零件i固定了一个切入点e_i。那么问题退化为一个标准的非对称旅行商问题（ATSP）吗？不，因为距离D(e_i, e_j)通常等于D(e_j, e_i)（欧氏距离对称），所以是对称TSP。我们的目标是找到访问{e_1, e_2, ..., e_N}这组点的一条最短哈密顿路径（如果起点固定，则是回路）。

如何求解这个TSP？对于N不太大（比如<50）的情况，我们可以使用模拟退火算法（SA）或遗传算法（GA）来获得一个优质解。这里我更推荐模拟退火，因为它实现简单，调整参数较少，适合作为核心优化器。

模拟退火算法设计要点：

状态：一个零件的排列序列π。
邻域操作：采用2-opt（交换两段路径）或随机交换两个零件的位置。2-opt在TSP中非常有效。
能量函数：即目标函数F(π) = D_start(e_{π1}) + Σ D(e_{πj}, e_{π(j+1)})。
降温计划：初始温度T0设置得足够高，使得初始的坏解也有较大概率被接受（例如，T0 = 1000 * F(initial)）。降温系数alpha通常取0.95到0.99。迭代次数（每个温度下的马尔可夫链长度）设为L = 100 * N。
终止条件：温度低于某个阈值T_min（如1e-6），或连续若干个温度下最优解未更新。

注意：这里存在一个常见的“坑”。我们计算的是点e_i到e_j的距离。但如果e_i和e_j不是简单的点，而是需要从轮廓点集合中动态选择呢？这就是我们下一步要解冻的。在第一步，我们强行固定了e_i，所以计算是直接的。这个固定策略为我们提供了一个基准序列。

3.2 策略二：固定访问序列，优化每个零件的切入点

现在，假设我们通过上一步得到了一个零件访问序列π。序列固定了，但每个零件的切入点e_{πj}还没定。这时，优化问题变成了一个动态规划（DP）问题，或者可以通过贪婪局部搜索高效求解。

问题描述：给定序列π1 -> π2 -> ... -> πN。对于零件πj，其切入点可以从其轮廓点集P_{πj}中任选一点p。我们需要为每个零件选择切入点，使得序列总空程最小。总空程S = D_start(p1) + D(p1, p2) + D(p2, p3) + ... + D(p_{N-1}, p_N)，其中pj ∈ P_{πj}。

这实际上是一个多阶段决策问题，非常适合用动态规划求解。

动态规划状态定义：

dp[j][k]：表示切割完前j个零件（即π1, ..., πj），并且第j个零件πj选择其第k个点作为切入点时，所累积的最小空程。
这里j从1到N，k是零件πj轮廓点的索引（0 到M_{πj}-1）。

状态转移方程：dp[j][k] = min_{t} { dp[j-1][t] + D( p_{π(j-1), t}, p_{πj, k} ) }其中p_{π(j-1), t}表示零件π(j-1)的第t个轮廓点。对于j=1（第一个零件）：dp[1][k] = D_start( p_{π1, k} )

最终答案：min_{k} dp[N][k]

通过DP，我们可以在O(N * M^2)的时间内找到给定序列下的最优切入点组合（M是平均轮廓点数）。如果M很大（如圆形离散化为360点），这个复杂度会很高。此时可以采用近似DP或贪婪法：对于零件πj，在固定π(j-1)切入点的情况下，选择πj轮廓上距离该点最近的点作为切入点。这种贪婪法一次遍历即可，复杂度O(N * M)，虽然不能保证全局最优，但效果通常很好，且可以作为DP的初始解或快速迭代工具。

3.3 策略三：迭代反馈与整体优化框架

现在我们把策略一和策略二结合起来，形成一个迭代优化框架：

初始化：
- 为每个零件i随机选择一个初始切入点e_i（例如，轮廓上的第一个点）。
- 或者，使用一个简单启发式：选择距离所有零件几何中心最近的点。
外层循环：迭代若干次（例如10-20次），或直到目标函数收敛。 a.阶段A：优化序列。固定当前所有零件的切入点{e_i}，将其视为TSP的城市，使用模拟退火算法求解最优访问序列π。 b.阶段B：优化切入点。固定上一步得到的最优序列π，使用动态规划或贪婪法，为序列中的每个零件重新计算最优切入点{e_i_new}。 c.更新与评估：用新切入点{e_i_new}替换旧的{e_i}。计算新配置下的总空程。如果优于历史最优，则更新历史最优解。 d.降温或扰动：模拟退火的外层也可以引入“温度”概念，以一定概率接受变差的切入点更新，避免陷入局部最优。或者，简单地在切入点更新后，对序列进行小幅扰动，重新进入阶段A。

这个框架将复杂的双层问题分解为两个相对简单的子问题，并让它们相互指导、迭代改进。在实际编程中，它非常有效，通常能在短时间内收敛到一个满意的解。

4. 代码实现详解与关键技巧

理论说完，我们来点实在的。以下是用Python实现上述核心框架的关键部分。我会用注释解释每一步，并分享一些调试和优化技巧。

4.1 数据结构定义

首先，定义零件（Part）类，用于存储轮廓信息和管理切入点。

import numpy as np import math import random import itertools class Part: def __init__(self, part_id, contour_points): """ 初始化一个零件。 :param part_id: 零件ID :param contour_points: 轮廓点列表，形状为 (n, 2) 的numpy数组，表示n个点的(x, y)坐标。 注意：点应该是闭合的，即 contour_points[0] 和 contour_points[-1] 应相同或非常接近。 """ self.id = part_id self.contour = np.array(contour_points) # 确保轮廓是闭合的 if not np.allclose(self.contour[0], self.contour[-1]): self.contour = np.vstack([self.contour, self.contour[0:1]]) self.num_points = len(self.contour) # 当前选择的切入点索引 self.entry_idx = 0 # 预计算轮廓周长（切割固定长度） self.perimeter = self._calculate_perimeter() def _calculate_perimeter(self): """计算轮廓周长。""" perimeter = 0.0 for i in range(self.num_points - 1): perimeter += np.linalg.norm(self.contour[i+1] - self.contour[i]) return perimeter def get_entry_point(self): """返回当前切入点的坐标。""" return self.contour[self.entry_idx] def set_entry_by_point(self, point): """给定一个坐标点，选择轮廓上离该点最近的点作为切入点。返回该点索引。""" distances = np.linalg.norm(self.contour - point, axis=1) self.entry_idx = np.argmin(distances) return self.entry_idx def get_closest_point_idx(self, point): """返回轮廓上离给定点最近的点的索引。""" distances = np.linalg.norm(self.contour - point, axis=1) return np.argmin(distances)

4.2 距离计算与目标函数

目标函数是空程，我们需要高效计算。

def euclidean_distance(p1, p2): """计算两点间欧氏距离。""" return np.linalg.norm(p1 - p2) def calculate_total_idle_distance(parts, sequence, start_point=np.array([0.0, 0.0])): """ 计算给定零件列表、访问序列和切入点选择下的总空程。 :param parts: Part对象列表 :param sequence: 零件索引的列表，表示访问顺序 :param start_point: 切割头起始位置 :return: 总空程 """ total_distance = 0.0 current_pos = start_point for part_idx in sequence: part = parts[part_idx] entry_point = part.get_entry_point() total_distance += euclidean_distance(current_pos, entry_point) current_pos = entry_point # 切割后，切割头仍在该点 return total_distance

4.3 模拟退火求解TSP（固定切入点）

这是策略一的核心。

def simulated_annealing_tsp(parts, initial_sequence, start_point, max_iter=5000, t0=100.0, alpha=0.95): """ 使用模拟退火求解TSP（零件访问序列）。 :param parts: Part对象列表，切入点已固定。 :param initial_sequence: 初始序列 :param start_point: 起始点 :param max_iter: 最大迭代次数 :param t0: 初始温度 :param alpha: 降温系数 :return: (best_sequence, best_distance) """ current_seq = initial_sequence[:] best_seq = current_seq[:] current_dist = calculate_total_idle_distance(parts, current_seq, start_point) best_dist = current_dist t = t0 n = len(parts) # 每个温度下的迭代次数，与问题规模相关 lk = n * 10 for iter in range(max_iter): for _ in range(lk): # 邻域操作：随机交换两个位置（另一种常用是2-opt，这里用简单交换） i, j = random.sample(range(n), 2) new_seq = current_seq[:] new_seq[i], new_seq[j] = new_seq[j], new_seq[i] new_dist = calculate_total_idle_distance(parts, new_seq, start_point) delta = new_dist - current_dist # 接受更差解的概率 if delta < 0 or random.random() < math.exp(-delta / t): current_seq, current_dist = new_seq, new_dist if current_dist < best_dist: best_seq, best_dist = current_seq[:], current_dist # 降温 t *= alpha if t < 1e-6: break return best_seq, best_dist

4.4 动态规划优化切入点（固定序列）

这是策略二的核心。这里实现贪婪法（最近点）作为示例，因为它简单高效。DP版本更精确但代码稍长。

def greedy_optimize_entry_points(parts, sequence, start_point): """ 贪婪法优化切入点：固定序列，每个零件选择离上一个点最近的轮廓点作为切入点。 :param parts: Part对象列表 :param sequence: 固定好的零件访问序列 :param start_point: 起始点 :return: 更新了切入点的parts列表，以及新的总空程 """ current_pos = start_point total_idle = 0.0 for part_idx in sequence: part = parts[part_idx] # 找到离current_pos最近的轮廓点，并设置为切入点 closest_idx = part.get_closest_point_idx(current_pos) part.entry_idx = closest_idx entry_point = part.get_entry_point() total_idle += euclidean_distance(current_pos, entry_point) current_pos = entry_point return parts, total_idle

4.5 主迭代框架

将以上模块组合起来。

def solve_cutting_path(parts, start_point=np.array([0.0, 0.0]), max_outer_iter=20): """ 主求解函数：迭代优化序列和切入点。 :param parts: Part对象列表 :param start_point: 起始点 :param max_outer_iter: 外层最大迭代次数 :return: (best_sequence, best_parts, best_total_idle) """ n = len(parts) # 初始化：随机序列，每个零件随机切入点 initial_sequence = list(range(n)) random.shuffle(initial_sequence) for part in parts: part.entry_idx = random.randint(0, part.num_points - 1) best_sequence = initial_sequence[:] best_parts = [part for part in parts] # 注意：这里需要深拷贝，简单起见用重新赋值切入点的方式 best_idle = calculate_total_idle_distance(parts, best_sequence, start_point) current_sequence = initial_sequence[:] current_parts = parts for outer_iter in range(max_outer_iter): print(f"Iteration {outer_iter+1}: Current best idle = {best_idle:.2f}") # 阶段A: 固定切入点，优化序列 # 注意：simulated_annealing_tsp内部计算距离时，使用的是parts当前的切入点 new_sequence, seq_idle = simulated_annealing_tsp( current_parts, current_sequence, start_point, max_iter=1000, t0=50.0 ) # 阶段B: 固定新序列，贪婪优化切入点 # 注意：greedy_optimize_entry_points会修改parts对象的entry_idx updated_parts, new_idle = greedy_optimize_entry_points( current_parts, new_sequence, start_point ) # 更新当前状态 current_sequence = new_sequence current_parts = updated_parts # 对象已更新 # 更新全局最优 if new_idle < best_idle: best_idle = new_idle best_sequence = new_sequence[:] # 保存最优状态下的切入点选择 for i, part in enumerate(current_parts): best_parts[i].entry_idx = part.entry_idx # 简单收敛判断：如果连续几次迭代没有改进，可以提前终止 # 这里省略，为了演示运行完整迭代 # 计算最终的总路径（空程 + 固定切割长度） total_cut_length = sum(part.perimeter for part in best_parts) total_path_length = best_idle + total_cut_length print(f"\nOptimization Finished.") print(f"Best sequence: {best_sequence}") print(f"Best idle distance: {best_idle:.2f}") print(f"Total cut length (fixed): {total_cut_length:.2f}") print(f"Total path length: {total_path_length:.2f}") return best_sequence, best_parts, best_idle

4.6 关键技巧与避坑指南

圆形离散化：如果零件包含圆形，务必将其离散化为足够多的点（如360）。点数太少会导致“最近点搜索”误差大，可能错过真正最优的切入点。在Part初始化时完成此操作。
起始点处理：切割头通常从一个“原点”或“换刀点”开始。我们的模型将D_start纳入目标函数是正确的。确保start_point参数设置正确。
模拟退火参数调优：t0（初始温度）和alpha（降温系数）需要根据问题规模调整。如果接受坏解的概率一开始就太低，算法容易陷入局部最优；如果降温太快，搜索可能不充分。一个实用的技巧是：让t0与初始解的目标函数值相关联（如t0 = 10 * initial_distance），并观察前几百次迭代中坏解的接受比例，将其调整到30%-50%左右为宜。
贪婪法与DP的选择：贪婪法（最近点）速度极快，O(N*M)，在迭代框架中作为默认选择。如果追求更高精度，可以在迭代的最后几轮或对最终序列使用动态规划（DP）进行精细优化。DP的O(N*M^2)复杂度在M较大时是负担，但可以尝试减少M（如对圆形只用72个点做DP优化）。
局部最优陷阱：我们的迭代框架本质上是交替优化，容易陷入局部最优。为了跳出，可以在外层循环中加入扰动机制：每隔几次迭代，随机改变几个零件的切入点，或者对序列进行一个较大的扰动（如随机反转一段序列），然后重新开始优化。
结果验证与可视化：一定要将最终路径画出来！使用Matplotlib将钢板、零件轮廓、空程（虚线）和切割路径（实线）可视化。肉眼观察路径是否交叉、是否明显绕远，这是发现模型或代码错误的最快方式。

5. 从模型到论文：解题思路的呈现与扩展

参加数学建模竞赛，光有代码和结果不够，还需要将你的思路清晰、逻辑严谨地呈现在论文中。针对这道题，论文的建模部分可以围绕以下几点展开：

问题重述与分析：强调问题的双层决策特性（序列与切入点），指出将其直接视为TSP的不足，引出分解与迭代的思想。
模型假设：明确列出切割头为点、空程定义、图形离散化等假设，使模型边界清晰。
符号说明：规范地定义文中使用的所有变量、符号，提升论文专业性。
模型建立：
- 整体模型：给出目标函数（总空程）的数学表达式。
- 子模型一（序列优化）：阐述在固定切入点下，问题如何转化为TSP，并说明采用模拟退火算法的理由（适用于组合优化、能处理中等规模问题、易实现）。
- 子模型二（切入点优化）：阐述在固定序列下，问题如何转化为动态规划或最近点选择问题，给出状态转移方程或贪婪策略。
- 迭代优化框架：用流程图展示“初始化 -> 序列优化 -> 切入点优化 -> 更新与判断”的完整流程，体现分治与迭代的思想。
算法步骤：用伪代码或清晰的步骤描述模拟退火、贪婪法/动态规划以及主循环的实现过程。
仿真结果与分析：
- 测试数据：自己构造几组不同数量、不同形状的零件数据（包括简单和复杂案例）。
- 结果展示：提供优化前后的空程对比数据表格。务必附上路径可视化图，这是最直观的证据。
- 灵敏度分析：探讨关键参数（如模拟退火的初始温度、降温系数，圆形离散化点数）对结果的影响。可以设计控制变量实验，用图表展示结果变化，体现研究的深度。
- 算法对比：可以简单对比纯贪婪算法（最近邻法）、只优化序列不优化切入点等方法，突出你模型的优越性。
模型评价与推广：
- 优点：分解思想降低复杂度，迭代框架保证解的质量，模型通用性强（可处理任意多边形和圆形）。
- 缺点：迭代法不能保证全局最优；对于零件数量极大（>500）的情况，计算时间可能较长。
- 推广：可扩展到考虑切割头加速度、不同切割速度、多切割头协同等实际场景。