简介：MATLAB实现粒子群算法（PSO）的完整代码包，面向计算机、电子信息、数学、物理、机械工程、土木工程等专业的大学生和研究生，适合毕业设计、课程设计或算法入门练习，以sum(x-0.5).^2为目标函数演示连续寻优过程，并绘制迭代曲线。资源共8个文件，包括6个M函数源文件、1个txt使用说明和1个docx程序说明，压缩包约22KB，结构紧凑、逻辑清晰。M函数采用模块化设计，将目标函数、初始种群生成、约束处理、速度限制、粒子解码与主流程分离开来，注释和参数说明详细，便于初学者修改目标函数、调整PSO参数，也可方便地迁移到其他优化任务；docx文档含程序说明和结果展示，txt文本介绍使用方法，能够帮助读者从代码和文档两个层面理解粒子群算法的实现细节。页面已有492人学习下载，适合作为PSO入门学习、实验教学和算法改进验证的参考实现。
1. 粒子群优化在 MATLAB 里的落地样板
当你面对一个不知道梯度、没有解析表达式的黑箱目标函数，粒子群算法往往比网格搜索快得多，也比遗传算法更容易调参。这个项目用不到 200 行 MATLAB 代码，把 PSO 的完整流程拆成了参数初始化、适应度计算、速度限幅、边界限幅和主循环几个独立文件，目标函数是简单的 y=sum(x-0.5).^2，迭代曲线一目了然。代码用的全是基础语法，从 MATLAB 2014a 到 2026b 的版本都能直接运行，不需要额外安装优化工具箱。对于做毕业设计、课程设计或刚开始研究群体智能算法的人来说，这套实现的价值不在那几句代码，而在“改目标函数、改参数、看效果”的路径非常短。
2. 粒子群算法原理与关键参数选型
2.1 位置-速度迭代模型与惯性权重
粒子群算法（Particle Swarm Optimization）的核心是让一群候选解在搜索空间中同步飞行。每个粒子拥有两个属性：位置 x 代表候选解，速度 v 代表搜索方向。每一代，粒子根据三部分信息合成新速度：上一代速度、自身历史最优位置 pbest、群体全局最优位置 gbest。用 MATLAB 向量化表达就是：
v = w * v + c1 * rand(size(x)) .* (pbest - x) + c2 * rand(size(x)) .* (gbest - x); x = x + v;第一项是惯性项，保留上一代飞行趋势；第二项是认知项，让粒子回到自己发现的好位置；第三项是社会项，让粒子向群体最优区域靠拢。三者合在一起，粒子既不会完全随波逐流，也不会只在自身周围打转。
w 的选择直接决定算法是“大范围探索”还是“小范围开采”。w 越接近 1，粒子越难改变既有方向，容易飞过最优点；w 越接近 0，粒子很快陷入当前最优附近的局部搜索，全局搜索能力差。常见做法是将 w 从 0.9 线性递减到 0.4，因为迭代初期粒子需要保持运动惯性覆盖全区域，后期则希望它能够稳定收敛到最优邻域。这套代码用的是w = w * damp的指数衰减，damp=0.99时 100 次迭代后 w 约为原值的 0.366，衰减速度比线性慢，适合需要长时间精细搜索的函数。
这里有一个经常被问到的点：pbest 是每个粒子自己的历史最优，gbest 是整个群体共享的。如果某个粒子一直没找到比 pbest 更好的位置，pbest 不变；但 gbest 一旦被其他粒子更新，所有粒子都会以新的 gbest 作为社会项目标，这就是粒子群能在群体层面出现“涌现”的原因。
2.2 学习因子、种群规模与边界约束
c1 和 c2 是认知项和社会项的加速系数。c1 过大会导致粒子频繁回到自身历史最优位置，群体内部沟通不足；c2 过大会让粒子过早聚到 gbest 周围，丧失探索能力。经典取法是 c1=c2=2，但在 Rastrigin 这类多峰函数上，我更倾向于用 c1 从 2.5 线性降到 0.5、c2 从 0.5 升到 2.5。原因是迭代后期需要粒子更信任群体信息，而不是个人经验。MATLAB 优化工具箱的 particleswarm 内部也有类似的自适应参数，但对外只暴露 SwarmSize、MaxIterations 等选项。如果你不想自己写主循环，也可以直接用工具箱验证结果，但自己做一遍的好处是可以实时观察粒子分布和收敛过程。
种群规模 nPop 的默认值不必随变量维度无限增加。对于 2 维到 10 维的问题，30 到 50 个粒子足够；维度超过 30 后再增加粒子数量效果不大，反而让每一步的适应度评估次数线性上升。对于 y=sum(x-0.5).^2 这个目标，甚至 10 个粒子都能在 50 代内找到大约 1e-10 的结果，不过为了曲线稳定，项目里取了 30 个粒子。
边界约束容易被新手忽略。速度更新公式没有“可飞范围”的限制，几次迭代后粒子位置就可能冲到 1e8 量级，之后的位置精度完全丧失。因此必须有 limspeedfun 和 limitfun 两道防线。速度上限一般取搜索区间宽度的 10%，比如区间 [-10,10] 则 vMax=2。边界处可以采用直接裁剪到边界，也可以把越界粒子位置重置为一个随机可行解。前者实现简单、收敛稳定，后者能增加多样性，但可能破坏已经形成的搜索方向，所以我通常选择裁剪。
下表是这套代码里主要参数的常用范围：
| 参数 | 常见范围 | 说明 |
|---|---|---|
| nPop | 20~50 | 粒子数，维度越高取越大 |
| maxIt | 50~300 | 最大迭代次数 |
| w | 0.9 衰减到 0.4 | 惯性权重 |
| c1 | 1.5~2.5 | 个体学习因子 |
| c2 | 1.5~2.5 | 群体学习因子 |
| vMax | (ub-lb)*0.1~0.2 | 速度限幅 |
上面这张参数表在后续的 main.m 里都能找到对应变量，改参数后不需要改动其他模块，这是模块化设计带来的最直接好处。
3. 模块化 PSO 核心代码：从目标函数到主循环
3.1 文件结构与职责划分
这套代码的文件划分很清晰，初学者拿到手之后不要急着跑，先对照下面的表确认每个文件在流程中的位置。
| 文件 | 职责 |
|---|---|
| myfun.m | 目标函数，返回 y=sum((x-0.5).^2) |
| genfun.m | 按给定维度生成初始粒子位置 |
| decodepsofun.m | 将粒子位置从编码空间映射到实际取值，常用于离散或混合变量 |
| limitfun.m | 对位置做边界约束 |
| limspeedfun.m | 对速度做最大最小值约束 |
| main.m | 主程序，完成参数设置、迭代循环、结果输出 |
这里要特别提一下 decodepsofun.m。很多初学者会忽略它：当变量是连续变量时，粒子位置本身就能直接代入目标函数，不需要解码。但如果你把 PSO 用在整数规划或二进制编码问题上，就需要单独写一个解码函数，把粒子的实数坐标映射到离散整数上去。这套代码里把解码独立成模块，意味着替换问题时你只需要改 myfun.m 和 decodepsofun.m，其他文件不用动。
3.2 主循环实现与全局最优更新
下面是从 main.m 提取并精简后的核心循环，可以直接贴到脚本里运行，配合前面的目标函数即可看到完整流程：
% main.m 核心循环 % 参数设置 nPop = 30; % 粒子数量 nVar = 2; % 变量维度 maxIt = 100; % 最大迭代次数 w = 0.9; % 初始惯性权重 damp = 0.99; % 每代衰减系数 c1 = 2; % 个体学习因子 c2 = 2; % 群体学习因子 % 初始化粒子结构体 particle = repmat(struct('pos', [], 'vel', [], 'cost', [], ... 'best', []), nPop, 1); gbest.pos = zeros(1, nVar); gbest.cost = inf; for i = 1:nPop particle(i).pos = genfun(nVar); % 生成初始位置 particle(i).vel = zeros(1, nVar); % 初始速度通常为零 particle(i).cost = myfun(particle(i).pos); % 计算适应度 particle(i).best = struct('pos', particle(i).pos, ... 'cost', particle(i).cost); % 个体最优 if particle(i).best.cost < gbest.cost gbest = particle(i).best; % 全局最优 end end % 迭代 for it = 1:maxIt for i = 1:nPop % 速度更新 vNew = w * particle(i).vel ... + c1 * rand(1, nVar) .* (particle(i).best.pos - particle(i).pos) ... + c2 * rand(1, nVar) .* (gbest.pos - particle(i).pos); particle(i).vel = limspeedfun(vNew); % 速度限幅 % 位置更新 pNew = particle(i).pos + particle(i).vel; particle(i).pos = limitfun(pNew); % 边界限幅 % 适应度计算与更新 particle(i).cost = myfun(particle(i).pos); if particle(i).cost < particle(i).best.cost particle(i).best.pos = particle(i).pos; particle(i).best.cost = particle(i).cost; if particle(i).best.cost < gbest.cost gbest = particle(i).best; end end end w = w * damp; % 惯性权重随迭代下降 history(it) = gbest.cost; % 记录每一代的最优值 end逻辑说明：每个粒子先用前一代速度、自身经验和群体经验三部分合成新速度，再限幅；位置更新后也必须经过边界限幅，然后计算适应度。只有个体最优改善时才会更新 pbest，同样只有当新 pbest 优于全局最优时才更新 gbest。gbest.cost初始化为 inf 是为了让第一代任意适应度都能覆盖它，这是一个容易被忽略的细节。
参数说明：rand(1,nVar)为每个维度独立生成随机系数，必须保持和位置向量相同的行向量形态。如果你不小心写成rand(nVar,1)，MATLAB 会做向量的隐式扩张或直接报错，最终结果精度会差几个数量级。粒子数量 nPop 和迭代次数 maxIt 决定总评估次数，值太小不收敛，值太大只是反复逼近同一个最优值，对精度提升有限。
3.3 速度限幅与边界限幅的实现
limspeedfun 和 limitfun 虽然短，却是 PSO 不在迭代中“飞散”的保险丝。我见过太多自己实现的 PSO 在 20 代后粒子位置全部变成 NaN，原因就是没有做速度限制。
% limspeedfun.m function v = limspeedfun(v) vMax = 1.0; % 速度上限，通常为变量范围的10% vMin = -vMax; v = min(max(v, vMin), vMax); end % limitfun.m function x = limitfun(x) lb = -10; % 下界 ub = 10; % 上界 x = min(max(x, lb), ub); end逻辑说明：限幅采用裁剪而不是随机重采样，是因为裁剪不会破坏粒子已经找到的有利方向，计算成本也最低。速度限幅的 vMax 如果取太大，限幅就失去意义；如果太小，粒子会像蜗牛一样爬，收敛极慢。边界限幅这里用的是“直接置为边界值”，还有一种策略是让粒子反射回可行域内部，但反射在多维空间里实现复杂，收益并不明显。
参数说明：lb 和 ub 应当与目标函数的真实定义域一致。对于 y=sum(x-0.5).^2 这个目标函数，理论上 x 的取值范围没有任何限制，但实际搜索必须给定一个有限区间，否则速度限幅的 vMax 无法确定。这个项目里默认把搜索区间设成 [-10, 10]，已经足以覆盖最优解 x=0.5 附近。
4. 迭代曲线绘制与收敛性验证
4.1 绘制适应度迭代曲线
第 3 章的主循环运行结束后，history 数组里存的就是每一代全局最优的适应度值。用下面的代码绘制收敛曲线：
% 绘制迭代曲线 figure; semilogy(1:maxIt, history, 'b-o', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('全局最优适应度值'); title('PSO 收敛曲线 (y=sum(x-0.5)^2)'); grid on;使用 semilogy 而不是 plot，是因为这个目标函数的适应度会从 10^2 量级快速掉到 10^-10 量级，线性坐标下后期曲线会变成一条贴地的直线，看不到波动。对数坐标能让你同时看到前期的大尺度下降和后期的微小改进。如果你还想同时画出 pbest 的均值或最差粒子，可以在主循环内添加一行avg_history(it) = mean([particle.cost]);，然后叠加绘制：
figure; semilogy(1:maxIt, history, 'b-o', 'LineWidth', 1.5); hold on; semilogy(1:maxIt, avg_history, 'r--', 'LineWidth', 1.2); xlabel('迭代次数'); ylabel('适应度'); legend('全局最优', '种群平均'); grid on;这样对比两条曲线的间距就能判断种群多样性。如果间距很小，说明粒子都聚在同一个点附近；如果间距一直很大，说明全局最优只属于少数粒子，可能陷入了局部极值。对于二维问题，还可以在最后一次迭代后用plot(gbest.pos(1), gbest.pos(2), 'rp')在目标函数等高线图上标出最优位置。
4.2 用不同参数对比验证收敛行为
验证 PSO 实现是否正常，最直接的办法是换参数跑几轮，对比收敛曲线。下面是一个可以自己复现的实验思路：分别用固定 w、线性递减 w、指数递减 w 各跑多次，统计最终的全局最优值量级。
| w 策略 | 典型收敛量级 | 说明 |
|---|---|---|
| 固定 w=0.7 | 1e-8 到 1e-10 | 后期速度无法减小，收敛慢 |
| 线性递减 0.9→0.4 | 1e-14 以下 | 推荐用于一般问题 |
| 指数衰减 damp=0.99 | 1e-14 以下 | 本项目采用 |
如果自己的实现与预期的量级差距超过 10 倍，先检查两个地方：一是随机种子是否固定，可以写rng(42);放在初始化前；二是速度更新时是否忘记使用.*逐元素乘法。向量方向不一致是 MATLAB 初学者最容易踩的坑，直接导致某些维度没有被正确搜索。
另外，熟悉 MATLAB 优化工具箱的话，可以用particleswarm验证自己的实现是否正确：
fun = @(x) sum((x - 0.5).^2); options = optimoptions('particleswarm', 'SwarmSize', 30, ... 'MaxIterations', 100, 'Display', 'final'); [x, fval] = particleswarm(fun, 2, -10, 10, options);这段代码里particleswarm的第三个和第四个参数是变量的下界和上界，返回的 x 是最优位置，fval 是对应适应度。如果自己的实现和工具箱结果一致，说明主循环逻辑正确。这个对照方法很适合用于课程设计报告中的结果验证部分。
5. 让 PSO 适配真实问题的三个技巧
真实问题的复杂性通常体现在三个方面：约束、多峰、混合变量。掌握下面三个技巧，就能从“复现算法”过渡到“改造成自己的优化工具”。
5.1 用罚函数处理约束
真实问题往往不是无约束的。最常见的是把不等式约束写成罚项加到 myfun.m 的返回值上：
function cost = myfun(x) y = sum((x - 0.5).^2); g = x(1) + x(2) - 1; % 示例约束 g(x) <= 0 if g > 0 cost = y + 1e6 * g^2; else cost = y; end end罚因子 1e6 要远大于目标函数的量级，否则约束被无视；但也不宜无穷大，否则适应度曲面会变成悬崖，粒子一越界就被弹回。通常罚因子取目标函数预期量级的 1000 到 10000 倍，需要根据实际输出调整。
5.2 按种群分布调惯性权重
指数衰减和线性衰减都属于开环策略，不看当前种群的分布。更有效率的做法是根据种群的“早熟程度”调整 w：当所有粒子的适应度方差很小时，说明都聚在一起，此时需要增大 w 跳出局部；当方差大时减小 w 加快收敛。计算方式：
f_avg = mean([particle.cost]); f_min = min([particle.cost]); sigma = std([particle.cost]); if sigma / (f_avg - f_min + eps) < 0.1 w = min(1.0, w * 1.2); else w = max(0.4, w * 0.98); end这个逻辑相当于给 PSO 加了一层反馈控制，在 Rastrigin、Griewank 这类多峰函数上比固定参数更有优势。
5.3 用解码函数处理离散变量
如果项目里包含离散变量，比如设备选型只能取 1、2、3，连续位置不能直接代入目标函数。这时 decodepsofun.m 就派上用场：把位置 x 四舍五入映射到整数集，然后再传给 myfun.m。注意速度更新仍然使用连续实数，只有解码后的值才进入目标函数，这样 PSO 的搜索动力不因量化而削弱。这也是这个项目把 decodepsofun.m 单独拆出来的最大价值。
最后补一个实操提示：跑完程序后，在命令窗口执行save('pso_result.mat', 'gbest', 'history')保存运行结果，方便写报告时复现曲线，不需要重新运行整个工程。
本文还有配套的精品资源，点击获取
