【数学建模学习笔记】启发式算法:粒子群算法

『AI先锋杯·14天征文挑战第5期』 10w+人浏览 445人参与

零基础小白看懂粒子群优化算法(PSO)

一、什么是粒子群优化算法?

简单说,粒子群优化算法(PSO)是一种模拟鸟群 / 鱼群觅食的智能算法。想象一群鸟在找食物:

  • 每只鸟(叫 “粒子”)不知道食物在哪,但能看到自己飞过的地方中 “最可能有食物” 的位置(个体经验);
  • 也能看到整个群体中 “最可能有食物” 的位置(群体经验);
  • 它们会根据这两个经验调整飞行方向和速度,慢慢靠近食物(最优解)。

PSO 就是用这种思路解决 “找最优解” 的问题,比如 “如何让某个函数的结果最小(或最大)”。

二、核心概念拆解

  1. 粒子:算法中的 “搜索者”,每个粒子有两个关键属性:

    • 位置:当前在搜索空间中的坐标(比如在三维空间中,位置可以用 (x1, x2, x3) 表示);
    • 速度:下一步移动的方向和快慢。
  2. 适应度函数:判断 “这个位置好不好” 的标准。比如我们想最小化函数f(x) = -10x1 - e^(-x2/10 - x3),那每个粒子的位置代入这个函数后,结果越小,说明位置越 “好”。

  3. 个体最优(p_best):每个粒子自己飞过的所有位置中,适应度最好的那个位置(“我之前飞过的地方里,这里最可能有食物”)。

  4. 全局最优(g_best):整个群体所有粒子飞过的位置中,适应度最好的那个位置(“大家飞过的地方里,这里最可能有食物”)。

  5. 位置和速度更新:粒子每次移动时,会根据以下规则调整速度和位置:

    • 速度 = 惯性(保持之前的速度) + 个体经验(向自己的 p_best 靠近) + 群体经验(向全局的 g_best 靠近);
    • 新位置 = 当前位置 + 新速度。

    就像鸟飞的时候,既不会完全忘记之前的方向(惯性),也会参考自己和同伴的最佳经验调整方向。

三、算法怎么跑起来?(用例子和代码说话)

假设我们要解决的问题是:最小化函数f(x) = -10x1 - e^(-x2/10 - x3),其中 x1 的范围是 0-1,x2 是 1-80,x3 是 0-120。

1. 初始化

先创建 100 个粒子(N=100),每个粒子有 3 个维度(x1, x2, x3,即 D=3);随机给每个粒子分配初始位置(在上述范围内)和初始速度(比如 - 1 到 1 之间)。

在代码中,初始化的实现如下:

# 导入必要的库
import numpy as np  # 用于数值计算
import random       # 用于生成随机数
import math         # 用于数学运算
import matplotlib.pyplot as plt  # 用于绘图

# 2. 初始化粒子的位置和速度
def initialize_particles(num_particles, num_dimensions, pos_low, pos_high, vel_low, vel_high):
    """
    参数说明:
    - num_particles:粒子数量
    - num_dimensions:变量维度(这里是3,因为有x1、x2、x3)
    - pos_low/pos_high:每个变量的位置范围(列表形式)
    - vel_low/vel_high:速度范围
    """
    # 初始化位置矩阵(行数=粒子数,列数=维度)
    positions = np.zeros((num_particles, num_dimensions))
    # 初始化速度矩阵
    velocities = np.zeros((num_particles, num_dimensions))
    
    for i in range(num_particles):  # 遍历每个粒子
        for j in range(num_dimensions):  # 遍历每个维度(变量)
            # 位置在[pos_low[j], pos_high[j]]范围内随机生成
            positions[i][j] = random.uniform(pos_low[j], pos_high[j])
            # 速度在[vel_low, vel_high]范围内随机生成
            velocities[i][j] = random.uniform(vel_low, vel_high)
    
    return positions, velocities

2. 定义适应度函数

适应度函数是判断位置好坏的标准,对于我们要最小化的函数,代码定义如下:

# 1. 定义适应度函数(目标函数)
# 这里要最小化的函数:f(x) = -10*x1 - e^(-x2/10 - x3)
def fitness(x):
    # x是一个列表,包含3个变量:x[0]是x1,x[1]是x2,x[2]是x3
    x1, x2, x3 = x[0], x[1], x[2]
    return -10 * x1 - math.exp(-x2/10 - x3)  # 计算函数值

3. 找初始最优

计算每个粒子的适应度(代入函数算结果);每个粒子的初始 “个体最优” 就是自己的初始位置;从所有粒子中挑出适应度最小的位置,作为 “全局最优”。

这部分在主函数中初始化个体最优和全局最优时实现:

# 初始化个体最优:一开始每个粒子的最优位置就是自己的初始位置
p_best_pos = positions.copy()  # 个体最优位置矩阵
# 计算每个粒子的初始适应度(目标函数值)
p_best_val = np.array([fitness(pos) for pos in positions])  # 个体最优值数组

# 初始化全局最优:从个体最优中找最好的
g_best_idx = np.argmin(p_best_val)  # 找到适应度最小的索引(因为我们要最小化函数)
g_best_pos = p_best_pos[g_best_idx].copy()  # 全局最优位置
g_best_val = p_best_val[g_best_idx]  # 全局最优值

4. 迭代优化(重复 100 次,即 M=100)

每个粒子根据 “惯性、个体最优、全局最优” 更新速度和位置;每次移动后,重新计算适应度,如果比自己之前的 “个体最优” 更好,就更新个体最优;再从所有个体最优中挑出最好的,更新全局最优。

速度和位置更新的代码实现:

# 3. 更新粒子的位置和速度
def update_particles(positions, velocities, p_best_pos, g_best_pos, w, c1, c2, pos_low, pos_high, vel_low, vel_high):
    """
    参数说明:
    - p_best_pos:每个粒子的个体最优位置
    - g_best_pos:全局最优位置
    - w:惯性权重(控制对之前速度的保留程度)
    - c1、c2:学习因子(控制个体/群体经验的影响程度)
    """
    num_particles, num_dimensions = positions.shape  # 获取粒子数和维度
    
    for i in range(num_particles):  # 遍历每个粒子
        for j in range(num_dimensions):  # 遍历每个维度
            # 生成0-1之间的随机数(增加搜索随机性)
            r1 = random.uniform(0, 1)
            r2 = random.uniform(0, 1)
            
            # 速度更新公式:惯性 + 个体经验 + 群体经验
            velocities[i][j] = (w * velocities[i][j]  # 惯性部分
                               + c1 * r1 * (p_best_pos[i][j] - positions[i][j])  # 向个体最优靠近
                               + c2 * r2 * (g_best_pos[j] - positions[i][j]))  # 向全局最优靠近
            
            # 限制速度在[vel_low, vel_high]范围内(避免速度过大)
            velocities[i][j] = max(min(velocities[i][j], vel_high), vel_low)
            
            # 位置更新公式:当前位置 + 新速度
            positions[i][j] += velocities[i][j]
            
            # 限制位置在[pos_low[j], pos_high[j]]范围内(避免超出变量定义域)
            positions[i][j] = max(min(positions[i][j], pos_high[j]), pos_low[j])
    
    return positions, velocities

5. 算法主函数

将上述步骤整合到主函数中,实现完整的 PSO 算法:

# 4. 粒子群优化算法主函数
def pso_algorithm(num_particles, num_dimensions, max_iter, pos_low, pos_high, vel_low, vel_high, w=0.9, c1=2, c2=2):
    """
    参数说明:
    - max_iter:最大迭代次数(搜索多少次)
    - 其他参数同上
    """
    # 初始化粒子位置和速度
    positions, velocities = initialize_particles(
        num_particles, num_dimensions, pos_low, pos_high, vel_low, vel_high
    )
    
    # 初始化个体最优:一开始每个粒子的最优位置就是自己的初始位置
    p_best_pos = positions.copy()  # 个体最优位置矩阵
    # 计算每个粒子的初始适应度(目标函数值)
    p_best_val = np.array([fitness(pos) for pos in positions])  # 个体最优值数组
    
    # 初始化全局最优:从个体最优中找最好的
    g_best_idx = np.argmin(p_best_val)  # 找到适应度最小的索引(因为我们要最小化函数)
    g_best_pos = p_best_pos[g_best_idx].copy()  # 全局最优位置
    g_best_val = p_best_val[g_best_idx]  # 全局最优值
    
    # 记录每次迭代的最优值(用于画图)
    best_values = [g_best_val]
    
    # 开始迭代搜索
    for iter in range(max_iter):
        # 更新粒子位置和速度
        positions, velocities = update_particles(
            positions, velocities, p_best_pos, g_best_pos, w, c1, c2, pos_low, pos_high, vel_low, vel_high
        )
        
        # 更新个体最优和全局最优
        for i in range(num_particles):
            current_val = fitness(positions[i])  # 当前位置的适应度
            # 如果当前位置比个体最优好(值更小),就更新个体最优
            if current_val < p_best_val[i]:
                p_best_pos[i] = positions[i].copy()
                p_best_val[i] = current_val
        
        # 从个体最优中找新的全局最优
        current_g_best_idx = np.argmin(p_best_val)
        current_g_best_val = p_best_val[current_g_best_idx]
        # 如果新的全局最优更好,就更新全局最优
        if current_g_best_val < g_best_val:
            g_best_pos = p_best_pos[current_g_best_idx].copy()
            g_best_val = current_g_best_val
        
        # 记录当前迭代的最优值
        best_values.append(g_best_val)
        # 打印迭代信息(每10次打印一次,避免输出太多)
        if (iter + 1) % 10 == 0:
            print(f"迭代第{iter+1}次,当前最优值:{g_best_val:.6f}")
    
    # 画收敛曲线(最优值随迭代次数的变化)
    plt.figure(figsize=(8, 5))
    plt.plot(best_values)
    plt.xlabel("迭代次数")
    plt.ylabel("最优目标函数值")
    plt.title("PSO算法收敛曲线")
    plt.grid(True)
    plt.show()
    
    return g_best_pos, g_best_val  # 返回最终找到的最优位置和最优值

6. 运行程序并查看结果

通过主程序调用 PSO 算法,设置参数并运行:

# 5. 主程序:运行PSO算法
if __name__ == "__main__":
    # 问题设置
    num_particles = 100  # 粒子数量(100个粒子一起搜索)
    num_dimensions = 3   # 变量维度(x1, x2, x3三个变量)
    max_iter = 100       # 最大迭代次数(搜索100次)
    
    # 变量范围:x1∈[0,1],x2∈[1,80],x3∈[0,120]
    pos_low = [0, 1, 0]
    pos_high = [1, 80, 120]
    
    # 速度范围(一般设为位置范围的10%-20%)
    vel_low = -1
    vel_high = 1
    
    # 运行PSO算法
    best_position, best_value = pso_algorithm(
        num_particles, num_dimensions, max_iter, pos_low, pos_high, vel_low, vel_high
    )
    
    # 输出结果
    print("\n优化结束!")
    print(f"找到的最优位置:x1={best_position[0]:.4f}, x2={best_position[1]:.4f}, x3={best_position[2]:.4f}")
    print(f"对应的最优目标函数值:{best_value:.6f}")

迭代结束后,全局最优位置就是算法找到的 “最佳解”。比如例子中最后得到的最佳位置可能是 (1, 1, 0),对应的函数结果约为 - 10.9048。

四、为什么要用 PSO?

  • 简单好懂:不需要复杂的数学推导,原理像 “群体找东西” 一样直观;
  • 灵活高效:能解决多变量、复杂函数的优化问题(比如最小化、最大化);
  • 容易实现:用 Python 等语言几行代码就能写出基本版本,如上述示例代码所示。

五、总结

粒子群优化算法(PSO)通过模拟群体中每个个体的 “学习”(参考自己和他人的经验)来搜索最优解,步骤简单、效果实用,是解决优化问题的常用工具。就像一群鸟通过互相 “借鉴” 飞行经验,最终找到食物最多的地方一样。

上述代码实现了粒子群优化算法,用于寻找目标函数f(x) = -10x1 - e^(-x2/10 - x3)的最小值(其中 x1、x2、x3 有明确的范围限制)。通过运行代码,可以直观看到粒子群算法如何通过群体协作逐步逼近最优解,适合零基础小白理解 PSO 的核心逻辑。如果想解决其他优化问题,只需修改 fitness 函数(目标函数)和 pos_low/pos_high(变量范围)即可。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值