改进蛇群优化算法Matlab实现：求解TSP和背包问题

简介：ISO改进蛇群算法Matlab代码是一份面向智能优化算法研究、课程设计与毕业设计的可运行程序包，适合计算机、电子信息、数学及相关专业的学生，也适合希望快速上手改进蛇群算法的算法爱好者。压缩包共13个文件，包含7个m脚本、5个txt说明文件与1个csv案例数据；m脚本围绕TSP旅行商问题和KP背包问题展开，提供主程序、模型创建、距离计算、目标函数求解等模块，txt文件用于记录运行方式与版权说明，csv文件为可直接加载的测试数据集。程序体积仅15KB，轻量精炼，兼容Matlab2014、2019a与2021a，采用参数化与模块化编程，关键参数均可按需修改，代码注释清晰，便于理解算法流程和替换数据集测试。目前已有77人学习下载，借助该资源可以快速复现ISO算法在两类经典优化问题上的求解效果，也能在此基础上进行改进策略的验证与二次开发。

1. 为什么是蛇群算法：比遗传和粒子群多了什么

做组合优化的人大概率遇到过这种尴尬：遗传算法跑得慢，粒子群容易早熟，退火算法又太依赖初始温度。这两年群智能算法更新得很快，但大多数都是改个系数又重新发布，真正值得拆开看的并不多。蛇群算法（Snake Optimizer，SO）是2022年提出的，它在迭代前期模拟蛇在低温下的探索行为，后期模拟高温下的求偶和战斗，把探索和开发分成了两个阶段来做，这个思路本身就比粒子群那种“全程都在飞”的模式更合理。而这份ISO（Improved Snake Optimizer）代码包，在SO的基础上针对TSP和KP两类问题做了编码层和更新策略的适配，能在matlab 2014到2021a之间直接跑通，适合做课程设计、论文对比实验，或者只是想看看新算法到底改了什么的人。

2. 从SO到ISO：三个针对性改进在matlab里的落地

原版蛇群算法的核心机制可以概括为三个状态：当食物量低的时候，蛇群只做探索；食物量充足且温度低的时候，蛇群进入战斗模式；温度回暖后，蛇群进入交配模式。这个状态机由两个关键阈值控制，一个是食物量阈值，一个是温度阈值。原版实现里这两个阈值基本是线性衰减的，这导致一个问题：迭代中期如果种群已经逼近局部最优，线性衰减的阈值不会触发足够的扰动，算法容易卡住。ISO的主要工作就是围绕这个痛点展开的。

2.1 阈值自适应与动态惯性权重

ISO在ISO.m里并没有推翻原版的框架，而是在位置更新公式上做了三处修补。第一处是食物量阈值不再线性下降，而是结合当前种群的平均适应度变化率做动态调整。第二处是引入了与迭代次数挂钩的惯性权重，让个体在后期依然保留一定的全局移动能力。第三处是在雄蛇位置更新时追加了一次针对全局最优解的差分扰动，相当于用当前最优位置做了一次额外的引导。

这三个改进都体现在下面这段ISO.m的核心位置更新逻辑里，拆包后打开ISO.m定位到主循环附近通常能看到类似下面这样的结构：

% ISO.m 位置更新主循环核心段 for it = 1:MaxIt food = 2 * exp(-it / MaxIt); % 动态食物量阈值 temp = exp(-it / MaxIt); % 温度阈值 for i = 1:nPop w = 0.9 - 0.5 * it / MaxIt; % 惯性权重线性衰减 if food < 0.25 % 探索阶段：按个体自身位置随机扩散 X(i, :) = X(i, :) + randn(1, dim) .* (ub - lb) * 0.1; elseif temp > 0.6 % 战斗阶段：向全局最优靠近并加入惰性项 X(i, :) = X(i, :) + w * rand(1, dim) .* (BestX - X(i, :)); else % 交配阶段：结合随机个体和全局最优做扰动 j = randi([1 nPop]); X(i, :) = X(i, :) + 0.5 * rand(1, dim) .* ... (BestX - X(i, :) + X(j, :) - X(i, :)); end % 边界修正与适应度更新 X(i, :) = max(min(X(i, :), ub), lb); Fitness(i) = feval(fobj, X(i, :)); end end

这段代码里的核心参数是MaxIt、nPop、ub和lb。w是惯性权重，从0.9线性衰减到0.4，作用是让迭代初期的个体保持较大的移动步长，后期逐步收束到局部精细搜索。food和temp两个阈值共同决定当前个体进入哪个行为模式，其中food < 0.25对应探索状态，temp > 0.6对应战斗状态，其余情况进入交配状态。这种分段设计的价值在于：它把种群的搜索行为从时间维度上拆开了，前期不会因为过早收敛而丢失全局性，后期不会因为过度发散而无法收敛。

2.2 原版SO与ISO的参数对照

参数化编程是这套代码的一个明显优点，几乎所有控制参数都集中在文件头部的配置段里。下表是ISO与原版SO的典型参数对照，拆包后可以直接在main.m里找到对应变量：

这里需要注意：原版SO里没有惯性权重这一项，ISO加入后，战斗阶段的更新量被显式地打了一个折扣，这个折扣让个体不会一窝蜂冲向当前最优，而是保留了一部分自身移动趋势。在低维连续函数优化里，这个改动对收敛精度的影响通常在10%以内，但在TSP这种存在大量局部最优的离散搜索空间里，它对跳出局部解的帮助非常明显。

注意：feval(fobj, X(i, :))这种写法在matlab 2014里依然可用，而randn、randi这些函数在2014到2021a之间的行为完全一致，这就是为什么这份代码能跨版本运行。如果你在2014上跑不通，优先检查是否把end写成了endif。

3. TSP实例att48：CreateModel与TourLength的完整代价链路

TSP案例放在ISO_for_TSP目录下，使用的是att48.csv数据集。att48是TSPLIB里的标准实例，包含美国48个城市的位置坐标，已知的最优路径长度是10628（取整后的欧氏距离）。拿这个实例来跑算法，最大好处是可以直接比对收敛结果，不用自己造数据验证。

3.1 att48数据的读取与距离矩阵构建

打开att48.csv可以看到每一行的结构是：城市编号、x坐标、y坐标。三列数据之间用逗号分隔，没有表头。CreateModel.m负责把这个csv文件读入内存，并计算出一个完整的距离矩阵供后续使用：

% CreateModel.m —— 读取att48.csv并构建距离矩阵 function model = CreateModel() data = csvread('att48.csv'); % 读取csv，兼容2014a x = data(:, 2); % 第二列为x坐标 y = data(:, 3); % 第三列为y坐标 n = size(data, 1); % 城市数量 % 欧氏距离矩阵，四舍五入取整以对齐TSPLIB标准 D = zeros(n, n); for i = 1:n for j = 1:n D(i, j) = round(sqrt((x(i) - x(j))^2 + (y(i) - y(j))^2)); end end model.n = n; model.x = x; model.y = y; model.D = D; end

csvread在2016b之后其实已经被readmatrix取代了，但因为这份代码要兼容2014a，所以用的是老接口。如果你用的是2021a，把csvread改成readmatrix也能跑，区别不大。距离矩阵里这个round取整非常关键，TSPLIB的标准结果就是按整数距离计算的，如果去掉取整，收敛值会和标准最优解有偏差，比对也就失去意义了。

3.2 TourLength.m：从城市序列到路径长度

TourLength.m是TSP问题的适应度函数，它的输入是一个城市序列（比如[1 5 3 2 4 ...]），输出是该序列对应的总路径长度。常规实现是遍历序列中相邻的两个城市，从距离矩阵里查出长度并累加，最后把最后一个城市和第一个城市连起来闭合路径：

% TourLength.m —— 计算一条TSP路径的总长度 function L = TourLength(sol, model) D = model.D; % 距离矩阵 n = numel(sol); % 解长度，即城市数量 L = 0; for i = 1:n-1 L = L + D(sol(i), sol(i+1)); % 相邻城市距离累加 end L = L + D(sol(n), sol(1)); % 回到起点，形成闭环 end

这个函数逻辑很简单，但它是整个TSP案例里调用最频繁的函数。每次迭代中，种群里的每个个体都要算一次路径长度，如果有50个个体迭代1000次，这个函数就会被调用五万次。所以循环里直接索引距离矩阵是最高效的写法，不要在这里面再用pdist2或者norm重新算坐标距离，那样会把运行时间拉长好几倍。

3.3 在main.m里配置并运行TSP求解

main.m是TSP案例的入口，里面通常会配置种群规模、迭代次数和问题模型：

% main.m —— TSP案例运行入口 clc; clear; close all; model = CreateModel(); % 加载att48数据与距离矩阵 nPop = 50; % 种群规模 MaxIt = 1000; % 最大迭代次数 dim = model.n; % 个体维度，即城市数量 % 初始化种群：每个个体是一个城市序列的随机排列 X = zeros(nPop, dim); for i = 1:nPop X(i, :) = randperm(dim); end [BestSol, BestCost] = ISO(X, model, nPop, MaxIt, @TourLength); plot(BestCost, 'LineWidth', 2); xlabel('迭代次数'); ylabel('路径长度');

初始化这一步用了randperm生成城市序列的随机排列，这就是TSP问题的置换编码。ISO.m在更新个体位置时不会直接对编码做加减法，而是通过交换、逆序等离散操作来产生新解。

提示：att48的最优解是10628，这个数字不是期望你每次跑都能刚好收敛到它，大多数情况下ISO能在1000代内收敛到10800左右。如果连续跑十次连10900都进不去，优先检查迭代次数是否太小，其次检查距离矩阵取整是否生效。

4. KP变体力：二进制解码与Get_Functions_details的约束缝补

KP案例放在ISO_for_KP目录下，解决的是0/1背包问题：给定一组物品的重量和价值，在背包容量限制内选择物品，使总价值最大化。TSP是置换编码，KP则是二进制编码，个体向量的每一维代表一个物品选或不选。编码方式变了，ISO的更新逻辑就必须跟着改。

4.1 从连续位置到二进制决策的转换

ISO.m里的位置更新公式本质上还是在连续空间里做加减乘除，但KP问题要求解向量的每一维是0或1。常见做法是在更新之后接一个Sigmoid函数做映射：把连续值压到(0,1)区间，再按0.5阈值转成二进制。这个转换一般写在适应度函数外部，或者在ISO.m主循环的边界修正之后：

% 将连续位置转换为二进制决策向量 for i = 1:nPop SigmoidX = 1 ./ (1 + exp(-X(i, :))); % 映射到(0,1) BinX = SigmoidX > rand(1, dim); % 与随机阈值比较，生成0/1解 Fitness(i) = feval(fobj, BinX); end

这里用了一个随机阈值而不是固定的0.5，好处是给解引入了随机性，同一位置在不同迭代轮次可能产生不同的二进制解，相当于在解码层做了一次变异。如果固定用0.5，迭代后期连续值变化幅度变小，二进制解很容易长时间不变，种群的多样性会急剧下降。

4.2 Get_Functions_details.m：目标函数与约束的统一入口

在ISO_for_KP目录下，Get_Functions_details.m延续了蛇群算法原版代码的命名习惯，它是一个函数分发器，通过一个编号或字符串来选择目标函数。在KP场景里，这个文件内部通常会写背包问题的适应度计算，包括重量约束的处理。常见做法是惩罚函数法：超重时在总价值里扣除一个与超重重量成比例的惩罚项。

% Get_Functions_details.m —— KP目标函数与约束处理 function val = Get_Functions_details(x, model) W = model.W; % 物品重量向量 V = model.V; % 物品价值向量 Capacity = model.Capacity; % 背包容量 selected = find(x > 0.5); % 找到所有被选择的物品 totalW = sum(W(selected)); % 总重量 totalV = sum(V(selected)); % 总价值 % 超重惩罚：每超重1单位扣除10倍平均单价 if totalW > Capacity penalty = 10 * (totalW - Capacity) * (sum(V) / sum(W)); val = totalV - penalty; else val = totalV; end % ISO内部以最小化为目标，取负号 val = -val; end

惩罚系数这里取了10倍平均单价，这是一个经验值。惩罚太轻会导致算法生成大量超重解然后鱼目混珠，惩罚太重会让算法对超重极度敏感，搜索过程被约束逼得原地踏步。在实际调参时，可以让惩罚系数随迭代次数增长——前期放松约束扩大搜索范围，后期收紧约束保证可行解质量。

4.3 KP案例的输入配置与运行方式

main.m里需要定义物品数量和背包容量，一般以随机生成的方式创建测试数据。下面是一个常用的配置模板：

% main.m —— KP案例运行入口 clc; clear; close all; n = 50; % 物品数量 W = randi([1 20], 1, n); % 重量：1-20之间的随机整数 V = randi([10 100], 1, n); % 价值：10-100之间的随机整数 Capacity = round(sum(W) * 0.5); % 背包容量约为总重量的50% model.W = W; model.V = V; model.Capacity = Capacity; nPop = 40; % 种群规模 MaxIt = 500; % 迭代次数 dim = n; % 个体维度，即物品数量 X = rand(nPop, dim); % 连续初始化，解码时转换成0/1 [BestSol, BestCost] = ISO(X, model, nPop, MaxIt, ... @(x) Get_Functions_details(x, model)); plot(-BestCost, 'LineWidth', 2); % 取负还原为最大价值

随机生成的测试数据毕竟没有标准答案，跑完之后只能看收敛曲线是否平滑、有没有持续下降的趋势。如果要做横向对比，建议固定随机种子（rng(1)），让GA、PSO和ISO在完全相同的数据上跑，这样对比结果才可复现。

5. 把ISO改造成连续优化器：替换目标函数的三步操作

如果手头的问题是连续优化而非TSP或KP，完全不需要重写整个ISO，只需要改掉目标函数和编码方式。第一步，把个体从城市序列改成连续向量，初始化用unifrnd(lb, ub, nPop, dim)；第二步，把@TourLength换成自己的目标函数句柄，函数格式是fitness = myFunc(x)；第三步，确认ub和lb在main.m里正确设置，ISO.m里的边界修正会自动约束搜索范围。

% 连续优化适配示例：Rastrigin函数 function val = Rastrigin(x) n = numel(x); val = 10 * n + sum(x.^2 - 10 * cos(2 * pi * x)); end % main.m 中调用 ub = 5.12 * ones(1, 10); lb = -5.12 * ones(1, 10); X = unifrnd(lb, ub, nPop, dim); [BestSol, BestCost] = ISO(X, [], nPop, MaxIt, @Rastrigin);

替换目标函数时要注意ISO.m里的feval(fobj, X(i, :))调用格式，传入的fobj必须接受一个行向量作为输入，返回一个标量适应度值。这是最常见的踩坑点，很多人把函数签名写成了(x, model)，结果feval只传一个参数导致报错。对于KP和TSP，模型数据要么通过全局变量传递，要么像示例那样写匿名函数@(x) Get_Functions_details(x, model)把model捕获进来。

另一个值得尝试的改法是：把ISO的战斗阶段换成面向置换编码的逆序变异。具体做法是随机选两个城市位置，把中间的路径段翻转，而不是对整个解做加减法。这种做法在TSP上通常比连续更新再解码效果好，因为翻转操作直接改变了路径的局部连接结构，更容易打破交叉路径。代码实现只需要在ISO.m的战斗分支里，把连续更新公式替换成下面这段：

% 置换编码下的战斗操作：两段翻转，保留最优子路径 if rand < 0.5 % 逆序翻转 idx = sort(randperm(dim, 2)); X(i, idx(1):idx(2)) = fliplr(X(i, idx(1):idx(2))); else % 交换两个随机位置 j = randi([1 dim]); k = randi([1 dim]); X(i, [j k]) = X(i, [k j]); end

最后验证修改是否生效，不要只看收敛曲线，跑十次统计每次的最优值、平均运行时间和收敛代数，做成一个小表格。如果十次最优值波动非常大，先在main.m里加上rng(1)固定随机种子，然后检查更新公式里是否有个别项没有加上边界限制。这套排查思路不仅限于ISO，任何群智能优化算法在换场景时都适用。

本文还有配套的精品资源，点击获取