鲸鱼算法去求解线性规划,听起来像是个有点奇怪的组合。线性规划是这个世界上最成熟、被研究得最透彻的优化问题之一,单纯形法、内点法加上各种商业求解器,已经把求解性能压榨到了极限;而鲸鱼优化算法(Whale Optimization Algorithm,WOA)是近几年才流行起来的群智能算法,靠模拟座头鲸捕食时的泡泡网策略,在解空间里到处"游动搜索"。这两者放在一起,很多人第一反应是:这图啥?
实话实说,如果只是为了求解一个标准的线性规划问题,直接调成熟的数学规划库是最优选择,完全没必要套一个启发式算法。但这个组合并不是没有价值。我在实际项目中遇到过太多"不标准"的优化问题:目标函数和约束条件里混进了非线性项,或者变量必须取整数,甚至干脆连模型本身都还在变化,这时候传统精确算法就不好使了。鲸鱼算法作为一个不依赖导数、不依赖凸性假设的通用搜索框架,可以直接顶上,而且它参数少、实现快、调起来也直观。
这个项目适合几类人:正在准备智能优化算法相关课程设计或者论文的同学,刚接触启发式算法、想知道它到底怎么和约束优化结合的人,以及有手头非线性或者混合整数优化需求、想找一个门槛低的基础求解框架的开发者。接下来我先把思路和原理拆明白,然后给出完整的Python实现和调试经验。
1. 为什么想到用鲸鱼算法去碰线性规划
1.1 线性规划的标准形式和精确算法的优势
线性规划的标准形式长这样:
min c^T x s.t. A x ≤ b, x ≥ 0其中 c 是目标函数系数向量,A 是约束矩阵,b 是约束右端项,x 是待求的决策变量。这类问题的核心性质是:可行域是一个凸多面体,最优解一定可以出现在某个顶点上。单纯形法沿着顶点一步步爬,内点法则从可行域内部逼近最优解,两者在凸性保证下都能高效收敛。
就拿我做过的某生产排产模拟项目X里的一个小模型为例,只有两个变量、两个约束,求解器在毫秒级别就能返回最优解,答案精度高到小数点后十位。这种场景下,智能优化算法确实完全不是对手。但问题在于,一旦模型变成非凸、离散、非线性,比如约束里出现 min、max、绝对值,甚至一个分段函数,顶点理论就失效了,精确算法要么没法建模,要么计算量爆炸。
这时候就需要一个"不管三七二十一先搜起来"的工具。鲸鱼算法就是干这个的,它不关心函数长什么样,只关心数值上能不能算出一个适应度分数。即使用它来解纯线性规划,也可以帮我们验证约束处理机制是否合理、种群搜索策略是否有效,以比较低的成本做算法验证。
1.2 座头鲸的捕食策略与优化搜索的对应关系
座头鲸在捕食磷虾群的时候,会下潜到鱼群下方,然后螺旋上升并吐出气泡形成气泡网,把鱼群逼向水面,最后在网中心张开嘴吞下猎物。鲸鱼优化算法把这个过程抽象成了三种行为:随机搜索、包围猎物、螺旋气泡网攻击。
三种行为正好对应优化算法里最基本的勘探与开发(Exploration and Exploitation)权衡。随机搜索负责大范围找新的区域,避免全军陷入一个局部最优;包围行为和螺旋更新负责在已知较好位置附近精细搜索,把解往上推。WOA最聪明的地方在于,它用一个收敛因子 a 和两个随机向量 A、C,在迭代过程中动态切换这三种行为,前期大范围勘探,后期精细开发。
我用一个生活类比解释:这就像是你在一个完全陌生的城市找一家最好吃的小店。随机搜索相当于每天随机坐公交去不同片区碰运气;包围相当于你听说某一片区有几家不错的店,于是围着那片区域挨家试;螺旋更新则像是在锁定的店铺附近细细品尝,从不同角度验证它是不是真的好。三种策略交替进行,既能保证不遗漏其他片区的黑马店,也能在已经发现的宝藏店周围深挖。
1.3 为什么我选了鲸鱼而不是粒子群或遗传算法
在这类项目里,可选的群智能算法不少,粒子群(PSO)、遗传算法(GA)、灰狼算法(GWO)都有人用。我最终选鲸鱼算法,有三个很实际的理由。
第一是参数少。粒子群要调惯性权重、个体学习因子、社会学习因子;遗传算法要处理编码方式、交叉率、变异率、选择策略。而WOA需要手工调的关键参数基本只有种群规模、最大迭代次数、螺旋形状常数 b,a 的衰减公式也是固定的。参数越少,调试成本越低,对新人越友好。
第二是更新逻辑统一。鲸鱼算法的位置更新只需要对当前解向量套用几个数学公式,没有交叉、变异这种离散操作,代码量非常少。对于线性规划这种连续且边界清晰的问题,"换个位置再评估一次"的思路特别顺畅。
第三是行为切换自然。收敛因子 a 从2线性递减到0,带动|A|从大于1逐渐变成小于1,算法自动从全局勘探过渡到局部开发,不需要人工干预。我在实验中测试过,同样的迭代次数下,WOA的收敛过程比PSO平滑很多,不会出现后期震荡剧烈的情况。
2. 鲸鱼算法核心机制拆解
2.1 从一条鲸鱼到一个候选解
在一个 n 维优化问题里,每条鲸鱼的位置就是一组 n 维坐标,表示一个候选解。我做的线性规划探秘项目中,决策变量是 x1 和 x2,那每条鲸鱼就是一个二维点 (x1, x2)。
初始化的时候,先读取每个变量的上下界,然后在边界范围内随机生成一批点。这就是初始种群。比如变量 x1 的取值范围是 [0, 10],x2 也是 [0, 10],那初始种群就是在这个正方形里随机撒几十个点,每个点代表一种生产方案的组合。
关键点在于:初始化阶段的随机性质量直接影响后续搜索效率。如果初始种群过于集中在一个角落,算法就要花更多迭代才能扩散到全局;如果太分散,又可能错过边界附近的可行域。我的做法是使用固定随机种子先跑一遍看分布,然后根据目标函数等值线的形状适当增加边界附近的初始点比例。
import numpy as np def init_population(pop_size, lb, ub): # lb和ub是形状为(dim,)的数组 dim = len(lb) pop = np.random.uniform(lb, ub, size=(pop_size, dim)) return pop2.2 包围、螺旋与随机的数学表达式
鲸鱼算法的核心是三套位置更新公式。假设当前最优位置是 X*,当前鲸鱼位置是 X,那么:
第一种是包围猎物:
D = |C * X* - X| X_new = X* - A * D其中 A = 2a * r1 - a,C = 2 * r2,r1 和 r2 是 [0,1] 区间的随机数。A 决定了这条鲸鱼向左逼近还是向右逼近,C 给 D 加了一层随机扰动。当 a 从2逐渐衰减到0时,A 的绝对值整体趋势越来越小,鲸鱼逐渐收缩到猎物周围。
第二种是螺旋气泡网攻击:
D' = |X* - X| X_new = D' * exp(b * l) * cos(2π * l) + X*l 是 [-1,1] 均匀分布的随机数,b 控制螺旋的疏密程度,通常取1。这个公式让鲸鱼在螺旋运动的同时逐渐逼近最优位置,模拟气泡网上升的收缩路径。
第三种是随机搜索:
D_rand = |C * X_rand - X| X_new = X_rand - A * D_randX_rand 是随机挑的一条鲸鱼,不是当前最优。这条公式的作用是当 A 绝对值较大时,让鲸鱼游向随机位置,保持全局勘探能力。
三种行为的选择逻辑是:生成一个随机数 p,如果 p < 0.5,走包围或者随机搜索路线,具体走哪条看 |A| 是否小于1;如果 p ≥ 0.5,走螺旋路线。这样设计的巧妙之处在于,螺旋更新和收缩包围是两种完全不同的搜索轨迹,同时保留两者可以增加种群行为的多样性,避免算法过早同质化。
2.3 勘探与开发如何随迭代自动切换
收敛因子 a 的衰减公式是:
a = 2 - 2 * t / Tt 是当前迭代次数,T 是最大迭代次数。刚开始时 a 接近2,随机向量 A = 2a * r1 - a 的绝对值大概率大于1,对应随机搜索行为,鲸鱼满空间乱跑。随着 a 越来越小,A 的绝对值逐渐小于1,算法开始收缩到最优解周边做精细搜索。
这个设计让我想到一个很实用的类比:打篮球时先退到中线附近观察整个球场,确定进攻方向后再运球进入禁区,最后在篮下不断调整出手角度。鲸鱼算法的前中后期就像一个球员逐渐从全局视野切换到篮下微调的过程,a 的线性衰减保证了切换不是瞬间发生的,给了算法充足的过渡时间。
实际项目里如果发现问题收敛太快或者太慢,我最先改的就是这个衰减公式。比如把线性衰减改成非线性衰减:
a = 2 * (1 - t / T) ** 1.5这个公式让前期衰减更快、后期衰减更慢,搜索重点更偏向后期的精细开发。遇到过目标函数有很多局部波动的情况,这么改效果立竿见影。
3. 实操:从零手写一个鲸鱼算法LP求解器
3.1 测试算例:一个两决策变量的线性规划
我用一个经典的两变量算例来演示。假设某工厂要安排两种产品的产量,目标是最大化利润:
max z = 3 * x1 + 2 * x2资源约束为:
x1 + x2 ≤ 4 2 * x1 + x2 ≤ 5 x1 ≥ 0, x2 ≥ 0用单纯形法很快能求出精确最优解是 x1 = 1,x2 = 3,最优利润 z = 9。我先用这个结果当"标准答案",再用鲸鱼算法去逼近它,看差距有多大。
这个算例虽然简单,但包含了两类典型约束:一个是主要资源的消耗约束,一个是组合约束,两个约束之间的交点正好就是最优解所在位置。用简单算例的好处是便于调试,如果算法连这个都跑不对,那一定是实现或者调参出了问题,而不是问题本身复杂度太高。
3.2 约束处理:罚函数这样写才靠谱
鲸鱼算法的搜索默认是无约束的,它只管让适应度值变小或变大,根本不认识什么 x1 + x2 ≤ 4。所以要把带约束的线性规划转化成算法能理解的形式,最常见的做法是罚函数法。
我的做法是构造一个新的适应度函数:
fitness = -z + penalty penalty = λ * (max(0, g1)^2 + max(0, g2)^2)其中 g1 = x1 + x2 - 4,g2 = 2*x1 + x2 - 5,约束违反多少,就加多少惩罚。注意负号是因为鲸鱼算法的核心逻辑是求最小值,而要最大化利润 z,所以取 -z。
罚函数系数 λ 的取值非常关键,我在这个项目里试过几组对比:
| λ 取值 | 现象 |
|---|---|
| 10 | 大量个体在可行域外游荡,结果偏小 |
| 100 | 勉强可行,但罚函数主导了搜索方向,收敛慢 |
| 1000 | 既能尽快回到可行域,又不会让目标梯度被惩罚项完全淹没 |
| 10000 | 惩罚过大,搜索变得很硬,容易卡在边界附近 |
最终我选了 λ = 1000。经验是:罚函数系数应该让"一个不可行解产生的惩罚"和"目标函数本身的量级"相当或者略大,这样搜索过程才会首先保证可行性,再去优化目标值。
在代码层面还需要一个小技巧:计算罚函数时,对不等式约束只惩罚超过 0 的部分,即 max(0, g_i),而不要对满足约束的部分进行任何评分,否则算法会过度偏向可行域内部,反而错过边界上的最优解。
3.3 完整代码与运行结果解读
下面给出完整的 Python 实现,里面保留了每个阶段的关键注释,方便直接照着改。
import numpy as np # 目标:max z = 3*x1 + 2*x2 # 约束:x1 + x2 <= 4, 2*x1 + x2 <= 5, x1>=0, x2>=0 def objective(x): return -(3 * x[0] + 2 * x[1]) # woa求最小,取负 def penalty(x, lam=1000.0): g1 = x[0] + x[1] - 4 g2 = 2 * x[0] + x[1] - 5 pf = lam * (max(0, g1) ** 2 + max(0, g2) ** 2) return pf def fitness(x): return objective(x) + penalty(x) def woa(pop_size=50, max_iter=200, b=1.0, lb=None, ub=None): if lb is None: lb = np.array([0.0, 0.0]) if ub is None: ub = np.array([10.0, 10.0]) dim = len(lb) # 初始化种群 pop = np.random.uniform(lb, ub, size=(pop_size, dim)) fitness_vals = np.array([fitness(ind) for ind in pop]) best_idx = np.argmin(fitness_vals) leader_pos = pop[best_idx].copy() leader_score = fitness_vals[best_idx] history = [] for t in range(max_iter): a = 2 - 2 * t / max_iter for i in range(pop_size): r1 = np.random.random() r2 = np.random.random() A = 2 * a * r1 - a C = 2 * r2 p = np.random.random() if p < 0.5: if abs(A) < 1: D = np.abs(C * leader_pos - pop[i]) new_pos = leader_pos - A * D else: rand_idx = np.random.randint(pop_size) X_rand = pop[rand_idx] D = np.abs(C * X_rand - pop[i]) new_pos = X_rand - A * D else: l = np.random.uniform(-1, 1) D_prime = np.abs(leader_pos - pop[i]) new_pos = D_prime * np.exp(b * l) * np.cos(2 * np.pi * l) + leader_pos # 边界处理 new_pos = np.clip(new_pos, lb, ub) pop[i] = new_pos fnew = fitness(new_pos) if fnew < fitness_vals[i]: fitness_vals[i] = fnew # 更新领导鲸鱼 best_idx = np.argmin(fitness_vals) if fitness_vals[best_idx] < leader_score: leader_score = fitness_vals[best_idx] leader_pos = pop[best_idx].copy() history.append(-leader_score) # 记录当前最优利润值 return leader_pos, -leader_score, history pos, best_val, hist = woa(pop_size=50, max_iter=200) print("最优解坐标:", pos) print("最优利润:", best_val)我实际跑出来的结果是 x 大约在 [1.0002, 2.9995] 附近,最优利润非常接近 9。从迭代历史来看,前 20 次迭代利润从很小一路爬升到 8 左右,20 到 60 次迭代慢慢逼近 9,后面基本稳定在 9,这个过程非常符合"前期勘探、后期开发"的预期。
如果当前环境中没有 NumPy,也可以直接用 Python 原生列表写一个简化版本,原理一样,只是循环慢一些。用 NumPy 的好处是批量计算方便,位置更新时能把整个种群向量化,后期扩展大规模算例时性能差距会很明显。
4. 调试记录:我在这个项目里踩过的四个坑
4.1 惩罚系数是命门,不是拍脑袋定的
第一次跑这个项目时,我偷懒把 λ 设成了 1,结果算法完全无视约束,最优解停在 x1 = 10, x2 = 0 这种边界点,利润看起来挺高,但资源约束被严重违反。后来我把 λ 提升到 10000,又发现算法像被"钉"在可行域边界上,每一步更新都小心翼翼,收敛速度慢得要命,迭代 500 次之后还没有逼近真实最优值。
教训是:调参的时候不要只看最终收益,还要监控两个指标。第一个是种群中"可行个体"的比例,第二个是当前最有个体对应的总约束违例量。如果违例量一直降不下来,加大 λ;如果违例量早就归零但利润进展很慢,适当调小 λ。我在代码里加了一步简单统计:
violate = max(0, pos[0] + pos[1] - 4) + max(0, 2 * pos[0] + pos[1] - 5)每迭代 10 次打印一次这个值,就能清楚地看到搜索是否在往可行域走。
4.2 迭代不收敛时先查这三件事
项目跑到后面,我故意把种群规模调成 5、迭代次数调成 30,结果最优解和精确值的相对误差飙到了 15% 以上。很多人这时候会开始怀疑算法不行,其实先别急着换算法,按顺序排查三件事。
一是种群规模。线性规划有两个决策变量,搜索空间是二维的,30 到 50 条鲸鱼就足够了;但如果是十个变量的实际问题,种群得放大到 100 到 200。种群太少,位置更新公式里的随机个体选择很容易把搜索引向死角。
二是最大迭代次数。可以用一个简单办法估算:看历史收敛曲线,如果利润还在明显上升就到了迭代末尾,那说明次数不够。比如上面示例代码里 60 次迭代已经接近最优值,200 次就非常稳妥。
三是随机种子。WOA 用到了大量随机数,不同种子跑出来的最优解会有小幅浮动。测试算法时固定随机种子,方便对比每次修改的效果;正式使用时再用多个种子跑多次,取最好结果。
4.3 边界反弹与死区的处理
鲸鱼算法理论上允许鲸鱼游出边界,但线性规划的决策变量有明确的上下界,有时候几个约束交汇的地方就落在某个边界上。我一开始没有做边界处理,导致不少鲸鱼被随机搜索带出可行区域,罚函数值巨大,浪费了大量迭代次数。
后来我在每次位置更新之后加上一行 np.clip(new_pos, lb, ub),把越界变量直接拉回边界刚内侧。这个操作看似简单,实际效果很显著。还有一个小技巧:在初始化时,我可以让 20% 的个体直接落在变量边界上,因为很多线性规划问题的最优解恰好在边界与边界的交点处,这么做等于给算法预先埋了"种子"。
4.4 别拿鲸鱼算法和单纯形法硬刚
这部分是我的真心话。如果问题是一个标准的、几百个变量的线性规划,闭着眼睛用成熟求解器就好,鲸鱼算法不管怎么调参,速度和精度都拼不过单纯形法或者内点法。鲸鱼算法的真正定位是"复杂问题的快速近似求解器"和"没有梯度信息时的兜底方案"。
| 场景 | 推荐算法 |
|---|---|
| 标准线性规划,变量连续,规模在几百以内 | 直接上 scipy.optimize.linprog |
| 带非线性目标或约束,无需精确解 | 鲸鱼算法先快速摸一遍 |
| 混合整数规划,规模不大 | 分支定界或者鲸鱼算法取整法 |
| 高维、大规模、非凸、无导数 | 鲸鱼算法 + 局部搜索精修 |
这不是算法的耻辱,而是工具的定位不同。就像家里有计算器,但还是有必要知道珠心算的原理,前者是生产工具,后者是理解和教学的工具。
5. 后续扩展:往前走一步更有意思
5.1 从标准LP到带非线性目标与约束的优化
把测试算例的目标函数从线性换成非线性,比如:
max f = 3 * x1 + 2 * x2 + sin(x1) * x2这个目标函数带有非线性项,传统单纯形法需要先把问题线性化,过程麻烦精度也受影响。鲸鱼算法几乎不用改代码,直接替换适应度函数就能跑,这也是我最初被它吸引的核心原因。
非线性约束也是一样,比如增加一个圆形约束:
(x1 - 1)^2 + (x2 - 1)^2 ≤ 9在罚函数里加上超出的量就行。这个圆形约束会让可行域变成一个非凸区域中间挖掉一块的形状,传统的顶点搜索彻底失效,而鲸鱼算法依然能通过搜索逐步逼近可行边界附近的最优解。我在测试中把过滤后的最优解画成散点图,能看到鲸鱼在可行域边界密集分布,这是算法在做边界开发的有力证据。
5.2 把鲸鱼算法改造成支持整数变量
很多实际排产问题要求决策变量是整数,比如生产台数、装箱数量。我做过一个小的模拟批次规划实验,处理方式是在位置更新之后加上一步:
new_pos = np.round(new_pos)但直接取整有个隐患:取整操作会破坏鲸鱼算法原本的连续运动轨迹,导致螺旋更新失去平滑性。我的改进方案是对整数变量使用罚函数软约束,不直接取整,而是在罚函数里加上:
integer_penalty = lam_int * (x - np.round(x)) ** 2这个惩罚项只会轻微地把解往整数点附近拉,不会像硬取整那样颠覆搜索路径。混合整数问题中一半变量连续一半变量整数时,这个方法尤其好用,连续变量正常更新,整数变量靠罚函数软性吸附到附近的整数值。
5.3 混合策略与并行计算
最后分享一个我在实际项目中用得很顺的技巧:把鲸鱼算法和成熟的局部搜索算法组合。具体做法是先用鲸鱼算法跑 100 次迭代,收敛到一个不错的解附近,然后将这个解作为初始点,交给 scipy 的序列最小二乘规划(SLSQP)做精细修正。前者的全局勘探能力弥补了 SLSQP 依赖初始点的缺点,后者的精确性弥补了鲸鱼算法收敛精度不足的问题,两者互补效果很好。
另外,鲸鱼算法每条鲸鱼的位置更新只依赖当前最优位置和随机个体,个体之间没有强耦合关系,非常适合并行评估。我用 multiprocessing 将种群分成若干组并行计算适应度,在 8 核机器上获得了接近 5 倍的加速比,处理 1000 维以上的高维优化问题时优势尤其明显。
这个项目整体做下来,我最大的体会是:智能优化算法和经典数学规划不是对立的,而是可以互相补充的两个工具箱。鲸鱼算法求解线性规划这件事,单看性能确实有点"杀鸡用牛刀",但它让我真正理解了约束怎么处理、勘探开发怎么平衡、参数为什么敏感这些问题。如果你也在做类似的算法探索,建议先从两变量小算例开始,把每步更新公式都验证一遍,再逐步扩大问题规模,最后你会发现,很多看似"不合理"的跨领域组合,反而藏着最有价值的思路。