PSO有约束优化：粒子群算法处理工程限制条件寻优实战

简介：解决带约束非线性优化问题的Matlab源码包，面向掌握一定优化算法基础的研究人员与工程师，采用Kennedy和Eberhart提出的粒子群算法，并借助Deb可行性规则处理等式/不等式约束，同时涵盖罚函数、约束投影等辅助策略。压缩包共8个文件：5个m脚本分别承担主程序、粒子群核心迭代、初始化、目标函数与绘图辅助功能；2篇PDF论文涉及液压缸优化设计和四个经典工程约束问题；另含1份Word文档，整体大小1.22MB，代码与文档配合紧密。目前已有741人学习下载。借助源码中的粒子位置更新、适应度计算、可行域校验等模块，再对照论文中的约束处理策略，可以快速理清全局最优解搜索流程，并理解不同可行度规则对收敛结果的影响。使用者可直接运行演示案例，也可依据自身工程场景调整参数与目标函数，降低从算法原理到实际应用的落地门槛。

1. PSO-有约束优化.zip 到底解决什么问题

很多工程师第一次拿到 PSO-有约束优化.zip 这个包时，会默认它是个普通优化器：扔一个目标函数进去，跑几百轮，出来一个最优值。但真正让这个包里代码和网上那些零散的粒子群 demo 拉开差距的，是它把“约束”当成了一等公民来处理——你做的刀具参数优化里主轴功率不能超，排产问题里交货期是硬性红线，这种带限制条件的极值问题才是工程常态。

这个包值得你花一个下午跑通，是因为它把粒子群搜索和约束处理组装成了一个可以直接替换目标函数的框架。你不用从零写粒子更新逻辑，只需要改自己的目标函数和约束函数。适合正在做参数寻优、资源调度、方案比选的工程师，也适合刚接触群体智能算法的学生用来对照课本理解收敛过程。一句话概括：它帮你解决的是“带限制条件的寻优”，而不是单纯的最小值搜索。下面从原理、复现步骤、调参和踩坑四个角度把它讲透。

2. 粒子群算法遇上约束：三种处理策略与选择理由

2.1 PSO 的核心机制与工程场景中的天然契合点

粒子群优化的灵感来自鸟群觅食。每只“鸟”代表候选解，在搜索空间里飞，飞行方向由两个经验支撑：自己的历史最佳位置 pbest 和整个群体的历史最佳位置 gbest。常见迭代公式是：

v = w * v + c1 * r1 * (pbest - x) + c2 * r2 * (gbest - x) x = x + v

四个关键部分各管一件事：w 是惯性权重，控制上一时刻速度保留多少；c1、c2 是加速系数，分别决定粒子向自身经验和群体经验靠拢的强度；r1、r2 是均匀随机数，给搜索加入探索性。这个公式没有梯度信息，只需要目标函数能给出输入对应的输出值，所以它能直接作用在“黑匣子”型目标函数上——比如外部仿真器的结果、试验测得的响应面，甚至查表得到的性能指标。

这类黑匣子场景在工程优化里非常常见。你没法写出加工时间关于切削参数的解析导数，但只要给一组参数能算出对应的加工成本和表面质量，粒子群就能跑。约束条件也一样，它不要求可导，只要求能判断“这个解是否越界、越了多少”。这个特性决定了 PSO 在有约束优化里特别吃香：约束函数的复杂度不影响算法迭代框架，只影响约束违反量的计算方式。

2.2 约束处理的三种主流策略

有约束优化问题的标准形式是：最小化 f(x)，同时满足不等式约束 g(x)≤0、等式约束 h(x)=0，以及变量上下界 lb≤x≤ub。处理约束的方式大体分三类：

这个包里最常见的设计是“罚函数为主、边界截断为辅”：变量越界直接拉回边界并把速度归零，工艺类不等式约束用外点罚函数折算进适应度。罚函数的标准写成：

fitness = f(x) + penalty_coef * max(0, g(x)) ** 2

这个形式的关键在于 max(0, g(x))：只有约束被违反时才产生罚项，平方是为了让违反轻微的点和严重违反的点拉开差距。等式约束 h(x)=0 则处理成 penalty_coef * h(x)^2，因为等式约束需要双边约束，绝对值在零点不可导，平方是通用选择。

选哪种策略，取决于约束本身的性质。如果约束允许一定程度的松弛（比如库存容量可以偶尔超一点，加急处理），罚函数法天生契合；如果约束是绝对安全红线（比如压力容器承压上限），可行性优先或严格罚函数更合适。我一般会先看约束个数和可行域占比：约束多且可行域窄时，纯罚函数容易把粒子全弹在边界外面，这时候混合策略更稳。

3. 跑通代码包：目录结构、最小命令与替换目标函数三步

3.1 压缩包内目录结构：先分清算法文件、测试函数和配置入口

解压之后不要急着运行，先花两分钟看清目录。典型布局如下：

PSO-有约束优化/ ├── src/ │ ├── pso.py # 粒子群算法主体 │ ├── fitness.py # 罚函数与适应度计算 │ └── constraints.py # 约束条件定义 ├── examples/ │ ├── run_demo.py # 最小运行示例 │ └── case_study.py # 带工业背景的案例 ├── config.yaml # 参数配置入口 └── README.md

用 tree 命令可以快速确认：

cd PSO-有约束优化 tree -L 2

如果没有 tree，Windows 下用dir /s，Linux 下用find . -maxdepth 2 -type f也能达到效果。确认这几个文件在不在，比直接运行更重要：pso.py 决定算法行为，fitness.py 决定罚函数怎么算，constraints.py 决定你的约束怎么写进框架。这三个文件是你要动手改的核心，其它都是辅助。

3.2 替换目标函数的最小步骤与验证流程

这类包的设计通常很统一：算法框架不动，用户填目标函数和约束函数。先用自带的 demo 验证环境没问题：

python examples/run_demo.py

跑通之后，按三步替换成你自己的问题。第一步，新建一个脚本 import 算法主体；第二步，定义目标函数和约束函数；第三步，初始化 PSO 对象并运行。示意如下：

import numpy as np from src.pso import PSO def my_objective(x): # 目标函数：例如两种原料的成本 return 3.0 * x[0] ** 2 + 2.0 * x[1] ** 2 def my_constraints(x): # 注意返回的是违反量，必须 >= 0 表示违反 return { "g1": x[0] / 5.0 - 1.0, # x0 不能超过 5 "g2": x[1] / 3.0 - 1.0, # x1 不能超过 3 } pso = PSO( pop_size=30, max_iter=200, w=0.7, c1=1.5, c2=1.5, bounds=[(0, 10), (0, 10)], penalty_coef=100, ) best = pso.run(my_objective, my_constraints) print(best.x, best.fitness)

这里最容易翻车的点是约束方向。不同作者写的包，约束函数返回值的含义可能相反——有些包约定了返回值大于 0 表示违反，有些则约定小于 0 表示违反。改之前先看 src/constraints.py 里对违反量的判断分支，或者在 demo 里故意设一个必违反的约束试一下输出。这个成本很低，但能省下后面排查的半小时。

替换完成后，强列建议用一个有解析解的小问题验证框架本身没被改坏。比如带线性约束的二次函数，理论上最优解位置可以直接手算出来。如果 PSO 跑出来的结果和手算值偏差超过 5%，说明你的目标函数或约束方向写错了，而不是算法不行。这一步是后面所有调参工作的基线。

3.3 用日志输出判断当前在“探索”还是“卡住”

很多包会打印每一代的当前最优值和约束违反量。如果没有日志，自己加三行记录：

for epoch in range(max_iter): pso.step() violation_sum = np.sum([abs(c) for c in pso.current_violations]) print(epoch, pso.gbest_fitness, violation_sum)

gbest_fitness 是到这一代为止的最优适应度，violation_sum 是当前种群所有粒子的约束违反量之和。这两条曲线的意义不同：fitness 下降说明在优化，violation_sum 下降说明在进入可行域。正常过程应该是两条曲线同时下降，最后 violation_sum 归零、fitness 收敛到稳定值。如果 fitness 一路降但 violation_sum 不为零，说明罚函数系数太低，粒子在用“违反约束换低目标值”；如果 violation_sum 归零但 fitness 长期不动，说明已经陷入局部最优或者罚函数系数过大导致搜索范围萎缩。

4. 调参才是收敛的关键：种群、惯性权重与罚函数系数的联动

4.1 五个核心参数的取值范围与物理含义

PSO 的参数不多，但相互耦合。核心参数表如下：

工程经验：w 从 0.9 线性递减到 0.4，是最稳妥的设置。前期高 w 让粒子大范围探索，后期低 w 让粒子在最有希望的区域精细搜索。c1 和 c2 通常取 1.5 左右的对称值，如果发现老是陷入局部最优，把 c1 调大到 1.8、c2 保持 1.5，增强每个粒子“走自己的路”的能力。

4.2 从默认参数到可行解的调试路径

不要一上来追求最优参数，先用一组保守默认值跑通：pop_size=40、max_iter=200、w=0.7 固定、c1=c2=1.5、penalty_coef=100。跑完之后看三件事：gbest_fitness 的下降曲线是否平滑、最终 violation_sum 是否归零、最终解的约束值离边界有多远。

如果 violation_sum 不为零，优先调大 penalty_coef。每乘 10 倍重跑一次，直到违反量消失。但要注意，罚系数不是越大越好：如果罚到 10000 以上粒子还在边界外徘徊，说明问题不在罚系数，而在初始种群——粒子生成时全落在不可行区域，罚函数又太陡，粒子根本没有机会看到可行域长什么样。这时候把初始种群改成“一半随机生成、一半按约束方向手动撒点”，效果立竿见影。

如果违规量已经清零但目标值不理想，问题多半在探索能力不足。把 w 提高到 0.9 或引入 w 线性递减策略，同时确认 max_iter 没设太少。一个需要记住的经验：PSO 的价值在于前 100 代快速逼近“好区域”，后 100 代才是在好区域里精细打磨。只跑 50 代，看到的基本都是半成品。

4.3 动态罚函数：后半程收紧约束的常见做法

固定罚系数的缺陷在于：前期约束严格会让粒子不敢探索，后期约束松弛又会让最终解落在边界之外一点点。更可靠的做法是让罚系数随迭代进度增长：

penalty_coef = 100 * (1 + epoch / max_iter * 99)

前期 penalty_coef 小，粒子可以相对自由地穿过不可行区域，找到目标值低的区域；后期 penalty_coef 接近峰值，约束变成硬性要求，把粒子压回可行域。这个策略在工程问题上几乎总是优于固定系数，因为它模拟了人做优化的习惯——先画大范围找方向，再逐步收紧条件。

动态罚函数另一个好处是减少了参数敏感性。固定的 penalty_coef 从 100 调到 1000 结果可能完全不同，但动态策略从“起止区间”这个层面控制，鲁棒性好很多。我一般把起止区间设成 10 到 1000 倍起步，如果最终解还有轻微违规，再把上限往上抬一个数量级。

5. PSO 有约束优化的常见坑与排查思路：现象、原因、解决

5.1 早熟收敛：所有粒子挤在一起却都不是可行解

现象：迭代到第 40 代左右 gbest_fitness 就不再变化，但 violation_sum 始终不归零，把最终解的约束值打印出来发现只是轻微越界。粒子位置在搜索空间里几乎重叠，多样性已消失。

原因：罚函数系数过大把非可行域变成了“高墙”，粒子在边界附近适应度差异极小，梯度信息被惩罚主导，失去了向可行域内部推进的方向。另一个常见诱因是初始种群全部落在同一片不可行区域，粒子根本没有见过可行解长什么样。

解决：把罚函数系数降回初始值的 1/10，让粒子能越过不可行区域，同时加大初始种群在搜索空间里的分布范围。如果种群已经挤成一团，别继续跑下去，重跑并换随机种子更省时间。也可以在迭代中期加一个重置机制：当 gbest 连续 20 代没更新，随机重置 20% 粒子的位置和速度。

5.2 罚函数系数失当：可行域边缘锯齿振荡与越界飞散

现象：迭代曲线是锯齿状，每一代 gbest 值上下跳动。有时候这一代找到一个可行解，下一代又跳到不可行区域里去。打印粒子位置时，发现大量坐标贴着边界值，比如 bounds 设的是(0,10)，有三分之一的粒子 x[0] 恰好是 10.0。

原因：罚函数系数太大导致所有非可行解的适应度都被压到无穷大附近，粒子一旦越过边界就被强行拉回，速度又没归零，于是来回震荡。越界粒子被拉回边界时速度还在，下一轮又冲出边界，形成锯齿。

解决：严格执行“越界截断”策略——坐标越界直接拉到边界，同时把该维度的速度置零，而不是只拉坐标不处理速度。代码上这样写：

for i in range(dim): if x[i] < lb[i]: x[i] = lb[i] v[i] = 0.0 elif x[i] > ub[i]: x[i] = ub[i] v[i] = 0.0

这段逻辑要放在位置更新之后、适应度计算之前。另外把罚函数系数改成动态增长，也可以缓解前期搜索被高墙弹回的问题。

5.3 速度越界与 NaN 传播的根因与修复

现象：某一次运行时，前面 30 代一切正常，第 31 代突然所有粒子的适应度都变成 nan，后面的日志全是 nan，程序也不报错但结果完全不可用。

原因：速度 v 没有做最大值限制，某次随机项叠加后粒子坐标飞到了极大值，比如 1e30 量级，目标函数内部如果算出指数或除法就直接溢出成 nan；或者目标函数自身在某组参数下抛出了除以零、对数负数的情况。nan 一旦出现会通过 pbest 和 gbest 传播到整个种群，且不可逆。

解决：强制限制速度范围。v_max 取搜索空间跨度的 10%~20%，例如变量范围是(0,10)，跨度 10，那 v 的每个分量限制在 [-2, 2] 比较稳妥。同时给适应度计算加一个保护分支：如果目标函数返回 nan 或 inf，直接把这个粒子的适应度置为 1e15，作为对照组垫底。另一个习惯是每次进入目标函数前检查输入是否在预设的“安全区间”内，越界就截断再进函数，从源头隔离掉除零和溢出。

5.4 黑匣子目标函数：单次仿真要数秒时的评估预算问题

现象：目标函数不是普通计算，而是一个外部仿真程序，跑一次要 5 秒。30 个粒子 × 200 代就是 6000 次仿真，等于 8 小时起步，根本不现实。

原因：没有把评估次数当成约束来规划。PSO 的评估次数 = pop_size × max_iter，这是这个算法最容易被忽略的“隐藏预算需求”。黑匣子场景必须牺牲种群规模或迭代次数。

解决：第一轮先用 15 个粒子 × 60 代跑通流程，确认目标函数与约束函数接口正确，输出不会闪退。然后调成 20 × 100 正式算一轮，记录最优结果。如果结果不理想，再加迭代数而不是加种群规模，因为 20 个粒子已经能维持基本多样性，而加迭代数不增加并行评估压力。常见的生产配置是 pop_size=20、max_iter=80，总评估 1600 次如果仿真时间 5 秒就是 2.2 小时，属于可接受范围。

6. 从测试用例到工程部署：约束松弛、离散变量与结果验证

6.1 把连续变量改成离散枚举的改法

工程问题里很多变量不是连续的，比如材料牌号只有 4 种、进给档位只有 6 档。PSO 本质在连续空间搜索，但可以做一个坐标映射层，把连续坐标映射到最近的离散档位。做法是：粒子在连续空间飞行，评估时把连续坐标按比例映射到离散档位：

levels = [1.0, 1.2, 1.5, 2.0] # 4 个离散档位 continuous_val = particle_pos[2] # 粒子在连续空间中的值 ratio = (continuous_val - levels[0]) / (levels[-1] - levels[0]) ratio = np.clip(ratio, 0.0, 1.0) idx = round(ratio * (len(levels) - 1)) discrete_val = levels[idx]

关键有两点：映射必须在每次评估前做，粒子实际看到的始终是离散值；粒子速度更新仍然基于连续坐标，不受映射影响。这个方案比直接在离散空间里改 PSO 公式简单得多，而且保持了连续空间搜索的平滑性。唯一要小心的是两个档位之间如果间隔很小，映射后会有大量粒子落到同一档位，种群多样性会下降，此时适当调大 w 能缓解。

6.2 用网格搜索校验 PSO 结果是否值得信赖

最后一步是验证 PSO 给出的解到底可不可信，不能看一条最优值就收工。标准做法是用完全独立的基准方法做交叉对比：对变量个数少的问题（2~5 维），网格搜索是够用的；变量多了就用随机搜索或者多起点局部搜索。

对比方法：用 3 个不同的随机种子跑 PSO，记录 3 次最优适应度和最优解坐标。把 3 次结果的差异算出来：如果 3 次最优值差异超过 10%，说明算法稳定性不足，需要调大种群或迭代；差异在 3% 以内说明稳定。再和网格搜索结果对比，PSO 的目标值应当优于网格搜索或至少相当。如果 PSO 反而比网格搜索差，多半是罚函数系数或约束方向出了问题，回去查而不是继续加迭代。

这套方法被我从测试函数带到生产环境的次数多了之后，我养成了一个习惯：不管约束看起来多简单，永远把约束函数单独写成一个文件。目的不是给代码看，而是留一个审计依据——出问题时你能说清楚目标函数和约束条件分别是什么。如果你也是第一次在工程上跑 PSO 有约束问题，建议先从这个包里的 demo 开始，改一个变量跑一遍，再逐步加约束。希望帮到你。

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