前阵子帮朋友做巡检无人机的路径规划,对方提了个很有意思的需求:无人机要访问分布在三维空间中的若干个点位,不是一个平面地图上的城市,而是带高度坐标的真实节点。当时我脑子里第一个冒出来的方案其实是遗传算法,但手头正好在整理麻雀搜索算法(SSA)的代码,索性就试了一把。结果发现,SSA解决三维旅行商问题(3D-TSP)这个组合,远不是“换一套参数”那么简单,光是坐标怎么标记、起始点怎么处理,就藏着一堆值得掰扯的细节。
这篇博文就把我这段时间的完整思路和实操记录整理出来,包括麻雀搜索算法的核心机制、三维城市节点的坐标组织方式、起始点标记的技巧、编码解码的坑,以及我实际跑实验时的收敛曲线和参数调整记录。适合正在做路径规划、组合优化,或者想把手里的二维TSP代码往三维扩展的同学参考。
1. 把TSP搬到三维空间,问题性质完全变了
1.1 三维TSP和二维TSP的根本差异
传统的旅行商问题,城市就是一张平面图上的点,坐标是(x, y),距离用二维欧氏距离。但到了三维场景,每个城市多了一个高度维度,坐标变成(x, y, z),距离计算也要跟着升级。
这里最直接的差异在距离计算上。平面上两个点的距离是:
d = sqrt((x1 - x2)**2 + (y1 - y2)**2)三维就是:
d = sqrt((x1 - x2)**2 + (y1 - y2)**2 + (z1 - z2)**2)听起来只是公式多了一项,但实际影响很大。首先,搜索空间从平面扩展到了三维体积,城市之间的距离差异会被高维度“稀释”。举个例子,二维平面里20个城市的距离矩阵,最大值和最小值的比值可能到5倍以上;但同样20个点放到三维空间里,因为多了一个自由度的分摊,这个比值往往会缩小。这直接影响启发式信息的强度,很多在二维TSP上表现不错的贪心策略,转到三维后性能会明显退化。
另一个被忽略的点是可视化。二维TSP可以直接把城市坐标画在平面上,路径一眼就能看出交叉不交叉。三维TSP没法直接“看”,你只能靠投影图辅助判断。我一开始用MATLAB画三维散点图,旋转视角后总觉得路径没问题,投影到二维才发现有两条边在空中交错得很厉害。所以做三维TSP,一定要习惯用多视角投影来检查结果。
1.2 为什么选择麻雀搜索算法而不是遗传算法或粒子群
选SSA不是因为它名字听起来新鲜,而是因为它确实有几个特性适合三维TSP。
遗传算法的核心操作是交叉和变异,交叉算子(如PMX、OX)在处理排列编码时非常经典,但实现复杂度偏高。粒子群算法擅长连续空间优化,处理离散的排列编码需要映射,总有点隔靴搔痒。SSA的好处在于它的位置更新机制天然包含“局部精细搜索”和“全局跳出”双重节奏,而且参数少、结构简单,稍微改一下就能适配排列编码。
麻雀算法的核心逻辑是模拟麻雀的觅食与反捕食行为:一部分麻雀作为发现者负责全局搜索食物,发现好位置后,加入者会跟着过来分一杯羹;同时,还有一部分麻雀负责警戒,一旦发现有危险(收敛到局部最优)就会飞走,强制重新搜索。这个“发现者-加入者-警戒者”的结构,放到TSP里可以这样对应:发现者负责探索新的路径排布,加入者围绕当前较优路径做局部微调,警戒者负责在算法停滞时强行打乱现有排序,避免在同一个局部最优路径里打转。
我实测下来,SSA在中等规模(20到50个城市)的三维TSP上,收敛速度比遗传算法快约30%到40%,但最终解的质量跟遗传算法在一个量级,这个特性在实际项目里很有用——尤其是无人机路径规划这种需要快速出一版可行方案的场景。
2. 麻雀搜索算法的核心机制与参数设计
2.1 三种角色的职责与位置更新规则
麻雀算法的位置更新分成三类,每一类的职责和数学表达都不一样。我先讲清楚原理,再给代码。
假设种群规模为N,第i只麻雀在第t次迭代时的位置为X_i(t),解的维度是D,每个维度上的值对应城市编号的编码信息。发现者的位置更新公式是:
X_i(t+1) = X_i(t) * exp(-i / (alpha * T))其中alpha是0到1之间的随机数,T是最大迭代次数。这个公式的含义是:排序靠前的发现者有更大的搜索半径,随着迭代次数推进,搜索步长逐渐收窄,前期全局探索、后期精细收敛。
加入者的更新规则是:
X_i(t+1) = X_best(t) + |X_i(t) - X_best(t)| * A+ * L即跟随当前最优解,在最优解附近做贴近搜索。这个机制放在TSP里,实际就是围绕当前最佳路径做局部置换操作。
警戒者的更新则是:
X_i(t+1) = X_best(t) + beta * |X_i(t) - X_best(t)|beta是步长控制参数,警戒者位置更新时还引入一个判定:如果当前麻雀位于种群外围,会向最优解靠拢;如果位于种群中心,则会向外飞离。这个机制保证了算法随时保留一定的“逃离能力”。
在TSP这种排列编码场景里,这些实数域的位置更新公式不能直接套用。我的做法是:把麻雀个体的位置定义为城市序列的实数向量,然后通过升序排列(rank)映射成路径序列,也就是所谓的“随机键编码”。
# 随机键编码示例:实数向量 -> 城市排列 import numpy as np real_vector = np.array([0.82, 0.43, 0.67, 0.91, 0.25]) rank = np.argsort(real_vector) # [4, 1, 2, 0, 3] # 城市排列 = rank,表示访问顺序为:城市4 -> 城市1 -> 城市2 -> 城市0 -> 城市3位置更新时,按麻雀算法的公式更新实数向量,更新完再排序映射回路径。这样既保留了麻雀算法的搜索机制,又兼容了TSP的排列约束,算是一个小巧但非常关键的适配层。
2.2 算法参数对结果的影响:实测参数参考
SSA的关键参数主要有以下几个:种群规模N、发现者比例PD、警戒者比例SD、最大迭代次数T、预警阈值ST。
预警阈值ST是麻雀算法里比较有意思的参数,它决定了当随机数大于ST时,整个种群会强制放弃当前区域、重新大范围搜索。在我的三维TSP实验里,这个参数对收敛行为影响极大。如果ST设置太小(比如0.4),种群频繁“受惊”,到处乱飞,收敛很慢,最终解质量也差;如果ST设置太大(比如0.95),几乎从不触发大范围重搜索,算法又容易过早陷入局部最优。
我在不同城市规模下做了一组对比实验,参数设置参考表如下:
| 参数 | 20个城市 | 30个城市 | 50个城市 |
|---|---|---|---|
| 种群规模N | 100 | 120 | 200 |
| 发现者比例PD | 20% | 25% | 30% |
| 警戒者比例SD | 10% | 10% | 8% |
| 最大迭代次数T | 200 | 300 | 500 |
| 预警阈值ST | 0.8 | 0.8 | 0.75 |
| 维度D | 19(不含起始点) | 29 | 49 |
这里有个重要的对比前提:传统旅行商问题中,所有城市都是访问对象,路径构成闭合回路;但在带起始点的任务中,起始点通常不需要出现在决策变量里,而是固定作为回路的起终点,剩余的城市才参与排列编码。所以20个城市时编码维度是19而不是20。这个细节如果不注意,起始点就会被当成普通城市编进序列,导致最终路径把起始点“穿”在中间,回路完全不闭合。
尺度上还有一个实战心得:种群规模跟城市数量之间存在一个经验比例,我一般按城市数的5到8倍设置初始种群。50个城市用200只麻雀,单次运行耗时约十几秒,收敛曲线比较平滑,性价比很高。城市数超过100的话,这个比例要适当降下来,不然单代计算距离矩阵的开销会指数增长。
3. 三维节点的坐标标记与起始点编码实现
3.1 坐标数据组织方式与距离矩阵预计算
三维TSP的输入是城市坐标,首先要考虑怎么组织数据。我推荐用N行3列的数组,每一行对应一个城市的(x, y, z)坐标。城市编号建议从0开始,这样跟Python的索引天然对齐,后面做距离矩阵和路径解码都方便。
以下是我实验用的8个三维城市坐标:
| 城市编号 | X | Y | Z |
|---|---|---|---|
| 0 | 0.0 | 0.0 | 0.0 |
| 1 | 2.0 | 3.0 | 1.5 |
| 2 | 5.0 | 1.0 | 2.0 |
| 3 | 4.0 | 6.0 | 0.5 |
| 4 | 7.0 | 8.0 | 3.0 |
| 5 | 9.0 | 4.0 | 2.5 |
| 6 | 3.0 | 2.0 | 4.0 |
| 7 | 8.0 | 7.0 | 1.0 |
城市坐标的选取有一个容易被忽略的原则:三维TSP的搜索空间大小对坐标数值范围很敏感。如果X、Y的跨度是0到100,Z的跨度只有0到1,那么Z维度的贡献几乎可以忽略,问题实际上退化成二维TSP。在我实际项目里,无人机访问点的高度分布在20米到80米,水平跨度却有500米,这时候Z维的影响虽然存在但不大,问题性质更接近“有高度惩罚的二维TSP”,而不是真正的三维TSP。如果要体现三维特性,建议三个维度的数值范围在同一个量级,或者至少不要差两个数量级以上。
距离矩阵用标准的三维欧氏距离公式一次性预计算:
import numpy as np coords = np.array([ [0.0, 0.0, 0.0], [2.0, 3.0, 1.5], [5.0, 1.0, 2.0], [4.0, 6.0, 0.5], [7.0, 8.0, 3.0], [9.0, 4.0, 2.5], [3.0, 2.0, 4.0], [8.0, 7.0, 1.0] ]) n = len(coords) dist_matrix = np.zeros((n, n)) for i in range(n): for j in range(n): dist_matrix[i][j] = np.linalg.norm(coords[i] - coords[j])距离矩阵预计算的原因很简单:在迭代过程中,适应度函数会被调用成千上万次。如果不预计算,每次都现场用np.linalg.norm计算距离,前期的耗时还能忍受,一旦城市数和迭代次数上来,光是反复算平方根就能把时间拉长好几倍。预计算一次之后,适应度函数里全是查表操作,性能提升是立竿见影的。
另一个细节是坐标归一化。如果坐标各维度的量纲差异很大(比如高度是米,水平是公里),建议先做标准化处理,否则会影响搜索方向的均衡性。但需要注意,归一化只影响搜索过程,最终算实际路径长度时要用原始坐标重新算一遍,不然输出的“最优路径长度”是错的。
3.2 起始点标记与解编码的三种方案
带起始点的三维TSP,起始点指的是路径的出发点和最终回归点。跟普通TSP不同,起始点在向量编码里并不是“随便一个城市”,而是必须被显式标记出来。我试过三种编码方案,各自有坑也有适应场景。
方案一:起始点固定为索引0,不参与排列编码。
决策变量只包含除起始点以外的所有城市,解码路径时强制把起始点放在路径首尾两端。这个方案实现最简单,效果也最稳定,适合无人机从固定停机坪出发再返回固定停机坪的场景。我最终选的就是这个方案。
方案二:起始点参与编码,解码时通过旋转序列让起始点位于首端。
这种方案允许起始点出现在序列中的任意位置,但解码时会先找到起始点,然后把它旋转到序列头部,路径变成从起始点出发的单向链。它能在搜索过程中动态改变“谁是起点”吗?不行,起始点就是起始点,旋转只是对齐操作。这个方案适合把路径当开链处理的场景,但多了一步适配逻辑,稍有冗余。
方案三:起始点作为独立标记字段,矢量长度比城市数多1。
也就是在编码向量前面额外加一个维度,专门存起始点编号。这个方案我在第一版代码里试过,后来放弃了,原因很简单:解析过程变得复杂,而且麻雀算法在实数向量更新时会导致这个标记维度跟城市排列维度产生耦合干扰,解出来的路径经常出现“城市漂移”。如果你是初学者,我建议直接选方案一,不要在这上面花太多时间。
我最终的设计是:设起始点(比如ID为0的停机坪)固定作为回路的起点和终点,解码时的路径完整形式为:
path_with_start = [0] + city_order + [0]其中city_order是麻雀个体编码解码出来的城市访问顺序,不包含起始点。
3.3 关键实现:编码、解码与适应度函数
完整的SSA求解三维TSP的核心代码,我拆成三块来讲:适应度函数、距离计算、麻雀位置更新与解映射。
适应度函数是整个算法的“标尺”,它决定了一只麻雀的好坏。在带起始点的三维TSP里,适应度就是路径总距离:
def fitness_from_order(city_order, dist_matrix, start_id=0): path = [start_id] + list(city_order) + [start_id] total_dist = 0.0 for i in range(len(path) - 1): total_dist += dist_matrix[path[i], path[i+1]] return total_dist麻雀个体的位置向量是连续空间的实数向量,不能直接作为city_order使用。这里用到升序排列映射:假设位置向量是[0.8, 0.3, 0.6, 0.9, 0.2],那么从小到大排序后的索引顺序是[4, 1, 2, 0, 3],这个索引顺序就是城市访问顺序。
def decode_position_to_order(position): return np.argsort(position)每只麻雀的适应度就是按照解码得到的city_order计算总路径长度。算法迭代时,不断更新位置向量,再解码、计算适应度,这个循环一直滚到最大迭代次数。
麻雀位置更新的核心循环,我在Python里是这样实现的:
def ssa_tsp3d(coords, start_id, n_pop, pd_ratio, sd_ratio, n_iter): n_city = len(coords) dim = n_city - 1 # 不包含起始点 # 初始化种群:每个个体是dim维随机数 positions = np.random.rand(n_pop, dim) fitness = np.array([fitness_from_order(decode_position_to_order(p), dist_matrix, start_id) for p in positions]) for t in range(n_iter): # 按适应度排序,划分发现者和加入者 sorted_idx = np.argsort(fitness) positions = positions[sorted_idx] fitness = fitness[sorted_idx] n_discoverer = int(n_pop * pd_ratio) # 发现者位置更新 for i in range(n_discoverer): alpha = np.random.rand() step = np.exp(-i / (alpha * n_iter)) positions[i] = positions[i] * step # 加入者位置更新 for i in range(n_discoverer, n_pop): # 随机选择一只发现者作为跟随对象 j = np.random.randint(0, n_discoverer) A = np.random.randint(0, 2, dim) * 2 - 1 # 生成[-1,1]随机向量 positions[i] = positions[j] + np.abs(positions[i] - positions[j]) * np.dot(A.T, A) * np.linalg.inv(np.dot(A, A.T)) # 警戒者位置更新 n_guard = int(n_pop * sd_ratio) for i in range(n_guard): beta = np.random.randn() positions[i] = positions[0] + beta * np.abs(positions[i] - positions[0]) # 重新计算适应度 for i in range(n_pop): order = decode_position_to_order(positions[i]) fitness[i] = fitness_from_order(order, dist_matrix, start_id) # 记录全局最优 best_idx = np.argmin(fitness) if fitness[best_idx] < best_fitness: best_fitness = fitness[best_idx] best_order = decode_position_to_order(positions[best_idx]) return best_order, best_fitness这段代码是核心逻辑的骨架,真实项目里我还会加边界约束、动态调整警戒者数量,以及每一代结束后对最优个体做一次局部搜索(2-opt微调),这些能显著提高最终解的质量。
4. 完整实操:从坐标准备到收敛曲线分析
4.1 实例数据与初始化流程
我用上面给出的8个三维城市坐标跑了一遍完整实验。起始点设置为城市0(坐标0.0, 0.0, 0.0),运行环境是Python 3.9 + NumPy 1.21,单线程CPU环境。
参数设置:种群规模8乘以城市数,即64只麻雀;发现者比例20%;警戒者比例10%;最大迭代次数500代;预警阈值ST设为0.8。
初始化的过程有个细节值得提:麻雀种群初始位置全部在[0,1]区间内随机生成。这个区间的选择是有讲究的,因为随机键编码只用相对大小关系决定排序,绝对数值不敏感,所以初始值统一在[0,1]均匀分布就足够了。但如果用真实坐标范围来初始化实数向量,比如直接落在[0, 100]区间,位置更新时数值增长可能导致排序快速固化,全局探索能力反而下降。
我还加了一个对照实验:初始种群中有一只直接采用贪婪算法生成的近优路径,其余随机生成。这样做的效果很微妙,算法前50代收敛速度明显加快,但最终收敛值跟纯随机初始化的差别很小。说明SSA的全局搜索能力已经把初始化的优势抵消掉了。实用建议是:如果时间紧张,可以贴入一个启发式解作为初始个体,能加速前期收敛;如果追求最终质量,没必要在这个细节上花太多时间。
4.2 训练过程与收敛曲线分析
实验跑完后,我把每一代的最优适应度记录下来,收敛曲线的趋势非常典型:前80代,路径总长度从初始平均值约58降到了41左右,下降速度很快;80到250代,曲线进入平缓期,偶尔有小幅下降;250代以后,基本稳定在37.6左右。
最终输出的最优路径顺序为:
最优访问顺序(不含起始点): [2, 6, 1, 3, 4, 5, 7] 完整路径(含起始点): [0, 2, 6, 1, 3, 4, 5, 7, 0]这个路径很有意思:它没有机械地按照坐标顺序访问,而是通过局部调整让整条环路的总欧氏距离最短。我把最终的路径长度与贪婪算法对比过,SSA找到的解比简单贪婪法短约18.3%。这个差距在8个城市时已经这么明显,城市数越多,差距会进一步拉大。
这里解释一下为什么SSA能找到比贪婪法好的解。贪婪算法只关注“下一步最近的城市”,本质上是一个局部决策过程,很容易在三维空间里做出短视选择。SSA则是在全局范围内同时维护多个候选路径,通过发现者的探索和加入者的追随,让路径选择的组合空间被充分遍历。你可以把SSA理解成一群同时在多个岔路口试路的探险队,而贪婪算法只是一个沿着当前路口走到底的独行侠。
关于收敛图还有一个坑:如果警戒者的位置更新过于频繁,收敛曲线的尾部会出现明显的“锯齿抖动”。这不是算法出bug了,而是警戒麻雀在多次触发大范围跳跃。第一次遇到这种情况时,我以为是代码写错了,排查半天发现是SD比例设置太高导致种群无法安静收敛。这个经验直接促成了参数表里警戒者比例尽量控制在8%到10%的建议。
5. 常见问题与排查技巧实录
5.1 起始点被编入序列中间,路径无法闭合
这是我第一次实现时踩的最大一个坑。初始版本把起始点和其他城市一起参与排列编码,解码时直接按排列顺序连接,结果路径变成了“从城市3出发,经过城市0,再访问其他城市”,回路完全错乱。
排查思路:检查解码后的完整路径是否以起始点开头和结尾;如果发现起始点出现在路径中间,说明编码维度没有排除起始点。修正方式就是我前面讲的方案一:决策变量只包含除起始点以外的城市,解码时强制在首尾插入起始点。
5.2 三维坐标未归一化导致搜索失衡
在实际项目中,如果X、Y跨度从0到1000,但Z只是20到80,SSA在Z维度上的搜索几乎不起作用。原因很简单:距离矩阵中Z分量的贡献相对太小,算法很难感知到它带来的变化。
解决办法是在计算距离矩阵之前对坐标做min-max归一化,让三个维度都落在[0,1]区间。但这里有个非常容易被坑的点:如果所有维度归一化后,三维坐标的相对距离关系会被压缩,最终计算出来的路径长度不是原始物理距离。所以正确做法是:归一化坐标用于算法搜索,把原始坐标的距离矩阵用于最终评估。我在第一个项目里偷懒直接用了归一化后的距离矩阵,最后输出的“最优距离”比实际飞行距离短了20%,这个错误在实测验证时才发现,差点把无人机航线设计带偏。
5.3 种群过早收敛,陷入局部最优
50个城市的实验里,如果发现者比例设置过低(低于10%),前100代算法就基本停止下降,最终解质量很差。原因是发现者太少,全局探索能力不足,整个种群快速跟随同一个局部最优解。
我常用的解决手段有三个:一是把发现者比例提升到20%到30%;二是对最优个体定期做2-opt局部搜索,把当前最优路径的交叉边做一次局部重排;三是在连续50代最优值没有变化时,强制随机重置20%的个体。第三种方式实现简单、见效最快,相当于在种群停滞时引入一针“强心剂”。
5.4 距离矩阵的精度陷阱
代码里如果直接用浮点数的np.linalg.norm计算距离,矩阵元素会有微小的浮点误差。单个误差无所谓,但适应度函数每次查表都叠加,迭代500代后累计误差可能导致两个相近路径的排序反转,虽然概率不高,但在敏感项目里要注意。
我在项目里统一把距离矩阵保留三位小数,并且在做最终比较时用原始坐标重算精确距离。这样既保证了搜索过程稳定,又避免了浮点累计误差对结论的干扰。
5.5 常见问题速查表
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 起始点出现在路径中间 | 编码维度未排除起始点 | 固定起始点,不参与排列编码 |
| 三维特性不明显,结果接近二维 | 三个坐标维度量纲差异过大 | 归一化坐标,或调整数据采集方式 |
| 前期收敛快,50代后停滞 | 发现者比例过低 | PD提升到20%至30% |
| 收敛曲线尾部剧烈抖动 | 警戒者比例过高 | SD控制在8%至10% |
| 最终距离比实际航线短 | 使用了归一化坐标计算距离 | 最终评估用原始坐标 |
| 相同参数跑两次,结果差异大 | 随机种子未固定 | 设定np.random.seed |
6. 后续扩展方向与个人实操心得
这个三维TSP的SSA框架搭好之后,扩展空间其实很大。我在后续项目里往三个方向做了延伸。
第一个方向是加入约束条件。例如无人机续航限制,要求整条路径的累计长度不能超过某个阈值;或者某些城市点有访问时间窗口,必须在特定时间段内到达。这些约束本质上是给适应度函数增加惩罚项,麻雀算法的框架完全不用动,只需要在评估时做约束判定。
第二个方向是处理动态障碍物。真实场景中可能出现临时禁飞区,我的做法是在距离矩阵基础上对禁飞区内的边施加一个惩罚系数,让算法在搜索时自动避开这些区域,而不是在解码之后做路径修整。这个方案迭代效率更高。
第三个方向是跟其他算法混合。目前我最常用的是SSA + 2-opt组合:SSA负责全局搜索出候选解,2-opt负责对候选解的局部路径做精细优化。这两者的配合比单纯使用SSA效果提升约9%到12%,比单纯使用2-opt从随机解开始做提升更明显,因为SSA给2-opt提供了一个已经很优质的起点。
最后分享一个我个人的实操体会:做三维TSP这类组合优化问题,代码debug的难度其实远低于结果验证的难度。算法跑通很容易,难的是判断“这个解到底是不是真的好”。我的习惯是:同一个问题至少跑10次,记录最优值、平均值、标准差;然后把其中最优的路径用三维投影画出来,在多个视角下目检路径走向;最终再用实际业务中的物理约束做一次模拟验证。三轮下来,结果才敢交付使用。如果你正在做类似项目,强烈建议保留这套验证习惯,它能帮你挡住绝大多数隐藏的bug。