数学建模_粒子群优化算法（PSO）_python

粒子群优化算法（PSO）笔记

1. 粒子群优化算法（PSO）概述

• 定义：粒子群优化算法是一种基于群体智能的优化算法，灵感来自自然界中鸟群或鱼群的觅食行为。PSO 算法通过在搜索空间中迭代地移动一组候选解（称为“粒子”）来搜索目标函数的最优解。
• 基本原理：每个粒子代表一个可能的解，并且通过迭代不断更新自己的速度和位置，最终收敛到全局最优解或局部最优解。

2. PSO 的关键概念

• 粒子（Particle）：每个粒子代表搜索空间中的一个位置，对应一个可能的解。
• 速度（Velocity）：粒子的移动速度，决定了粒子在搜索空间中的移动方向和距离。
• 位置（Position）：粒子在搜索空间中的当前坐标，对应于当前解。
• 个体最佳位置（pbest）：粒子自身找到的历史最佳位置。
• 全局最佳位置（gbest）：当前所有粒子中找到的最佳位置。

3. PSO 更新公式

• 速度更新公式 ：

• 位置更新公式：
  

4. Python 中的 pso 函数

• 来源：pso 函数来自于 pyswarm 库，用于实现 PSO 算法的优化过程。
• 功能：优化用户定义的目标函数，返回最优解的位置（best_position）和最优目标函数值（best_cost）。

5. pso 函数的关键参数

• Objective_function：目标函数，PSO 算法试图最小化或最大化的函数。
• lower_bound 和 upper_bound：搜索空间的上下界，定义粒子位置的取值范围。
• swarmsize：粒子群的数量，决定了算法同时探索多少个解。
• maxiter：最大迭代次数，控制算法执行的最大步数。
• debug：调试模式，设置为 True 时会输出详细的迭代信息，帮助理解优化过程。

6. 结果解释

• best_position：粒子群优化算法找到的最优解的位置（坐标）。如果是多维优化问题，best_position 是一个数组，包含每个维度的最优解。
• best_cost：目标函数在 best_position 处的值，通常是最小化或最大化的目标函数值。

7. PSO 的优势

• 全局搜索能力强：PSO 能够在较大的搜索空间内进行全局搜索，避免陷入局部最优解。
• 实现简单：PSO 算法结构简单，易于实现且计算效率较高。
• 适应性强：PSO 适用于多种优化问题，包括连续型和离散型优化问题。

8. PSO 在 Python 中的应用示例

best_position, best_cost = pso(Objective_function, lower_bound, upper_bound, 
                               swarmsize=n_particles, maxiter=max_iter, debug=True)

• 用途：通过调用 pso 函数，对目标函数进行优化，找到最优解的位置和对应的目标函数值。
• 典型场景：参数优化、函数最优化、模型调参等问题中广泛应用。

总结

PSO 是一种强大的优化算法，通过模拟群体智能来寻找最优解。它实现简单、全局搜索能力强，适合应用于各种优化问题。通过合理设置粒子数量、迭代次数等参数，PSO 可以在不同的应用场景中找到接近最优的解。

代码示例：

import numpy as np
import matplotlib.pyplot as plt
from pyswarm import pso

list_l = ['花叶类', '花菜类', '水生根茎类', '茄类', '辣椒类', '食用菌']

predit_buy = [
    [3.285864, 3.2921748, 3.2889733, 3.285188, 3.2851105, 3.2964268, 3.2876368],
    [7.7414317, 7.763459, 7.814592, 7.794937, 7.747068, 7.810813, 7.7633805],
    [12.018651, 11.912668, 12.027704, 11.941088, 11.92054, 12.118359, 11.972251],
    [4.5562034, 4.601929, 4.5483465, 4.549116, 4.532483, 4.539543, 4.601603],
    [3.7067149, 3.65774, 3.6644902, 3.6755412, 3.6658049, 3.6834998, 36471841],
    [4.015016, 4.075036, 4.025253, 4.0783653, 4.0397897, 4.0211616, 4.0569587]
]

predit_sale = [
    [190.75572, 189.94437, 189.2342, 190.34938, 189.08669, 188.56415, 190.10588],
    [28.618061, 28.872581, 28.873682, 28.74203, 28.776909, 28.86964, 28.661997],
    [26.962397, 27.805391, 27.65219, 27.210875, 27.88252, 27.764929, 27.24346],
    [32.29913, 31.795496, 32.425, 31.649815, 31.683603, 31.381622, 31.967655],
    [102.14867, 101.67641, 102.01936, 102.55891, 102.07538, 102.73362, 102.2766],
    [56.307552, 57.39569, 58.130955, 56.80816, 56.4629, 57.186737, 57.953686]
]

predit_omega = [0.7]*6  # 折扣
predit_gama = [0.1283, 0.1551, 0.1365, 0.0668, 0.0924, 0.0945]  # 损耗率

day = 6

ini_pos = [0.6, 0.4, 0.3, 0.6, 0.9, 0.6, 150, 30, 20, 20, 100, 60]
ini_pos = np.array(ini_pos)
n_particles = 1000
n_dimensions = 12
lower_bound = np.array([0.3, 0.2, 0.2, 0.3, 0.3, 0.3, 10, 0, 0, 10, 10, 20])
upper_bound = np.array([1.2, 0.98, 0.8, 0.9, 1.2, 1.0, 450, 90, 75, 60, 300, 250])

# 增加迭代次数
max_iter = 300

# 初始化成本历史记录
cost_history = []

def modify(y, x, idx):
    if idx == 0:
        sale = -1.91 * x + 1.03 * y + 3.82
    elif idx == 1:
        sale = -0.25 * x + 0.94 * y + 1.07
    elif idx == 2:
        sale = -0.03 * x + 0.97 * y - 0.98
    elif idx == 3:
        sale = -0.07 * x + 0.99 * y + 1.70
    elif idx == 4:
        sale = -0.029 * x + 0.99 * y + 2.54
    elif idx == 5:
        sale = -0.51 * x + 1.01 * y + 2.32
    return max(0, sale)

def Objective_function(x):
    profit = 0
    for idx in range(6):
        alpha = x[idx]  # 利润率
        beta = x[idx + 6]  # 订货量
        buy = predit_buy[idx][day]  # 预测购买价格
        sale_lstm = predit_sale[idx][day]  # 预测销售量
        omega = predit_omega[idx]  # 折扣系数
        gama = predit_gama[idx]  # 损耗率
        sale_price_normal = buy * (1 + alpha)  # 正常销售价格
        sale_price_discount = sale_price_normal * omega  # 折扣销售价格
        sale_modify = modify(sale_lstm, sale_price_normal, idx)  # 根据价格调整后的销售量
        w1 = sale_modify * (1 - gama) * sale_price_normal + sale_modify * gama * sale_price_discount
        w2 = beta * buy  # 总成本
        if beta <= sale_modify:
            profit += (w1 - w2) - (sale_modify - beta) * 20  # 销售损失惩罚
        else:
            profit += (w1 - w2)  # 利润计算
    cost = -profit  # 最小化负利润（最大化利润）
    cost_history.append(cost)
    return cost

# 使用PSO优化
best_position, best_cost = pso(Objective_function, lower_bound, upper_bound, swarmsize=n_particles, maxiter=max_iter, debug=True)

# 设置全局样式
plt.style.use('seaborn')

# 增加 DPI 以提高图像分辨率
fig, ax = plt.subplots(figsize=(8, 6), dpi=300)

# 绘制目标函数（-利润）变化曲线
ax.plot(cost_history, label="Objective Function (-Profit) Change", color='b', linewidth=2)

# 设置网格线
ax.grid(True)

# 设置标题、标签和图例
ax.set_title("Particle Swarm Optimization Objective Function Curve", fontsize=14)
ax.set_xlabel("Iteration", fontsize=12)
ax.set_ylabel("Objective Function Value", fontsize=12)
ax.legend(loc="upper right", fontsize=10)

# 调整坐标轴刻度和范围
ax.tick_params(axis='both', which='major', labelsize=10)

# 添加LaTeX形式的标签 (可选)
ax.set_xlabel(r'Iteration Number', fontsize=12)
ax.set_ylabel(r'Objective Function Value ($10^7$)', fontsize=12)

# 标注收敛点或其他关键点 (可选)
final_cost = cost_history[-1]
ax.annotate(f'Final Cost: {final_cost:.2e}', xy=(len(cost_history) - 1, final_cost),
            xytext=(len(cost_history) * 0.7, final_cost * 1.5),
            arrowprops=dict(facecolor='black', shrink=0.05),
            fontsize=10, color='red')

# 导出高质量图片
plt.savefig('pso_optimization_curve.png', bbox_inches='tight')

# 显示图表
plt.show()

def pre_sale_mount(x1, x2):
    list_sale = []
    for idx in range(6):
        x = x1[idx]
        y = x2[idx]
        sale = modify(y, x, idx)
        list_sale.append(sale)
    return list_sale

def count_how_much(x, y):
    list_how_much = []
    for i in range(6):
        list_how_much.append((x[i] + 1) * y[i])
    return list_how_much

print("最大利润：", -best_cost)
print("进货量（变量）：", best_position[-6:].tolist())
print("\n预测进价：", [row[day] for row in predit_buy])
sale_price = count_how_much(best_position[:6].tolist(), [row[day] for row in predit_buy])
print("销售价格：", sale_price)
print("利润率（变量）：", best_position[:6].tolist())
print("修正之后的销量：", pre_sale_mount(sale_price, [row[day] for row in predit_sale]))