每日算法 — 使用java实现粒子群优化：群体智能启发式全局搜索算法

可莉爸爸 AI算法

引言

粒子群优化（Particle Swarm Optimization，PSO）是由Kennedy和Eberhart于1995年提出的一种基于群体智能的启发式优化算法。灵感来源于鸟群觅食和鱼群游动的集体行为，PSO通过模拟群体中个体之间的协作与信息共享，在复杂搜索空间中寻找全局最优解。与遗传算法相比，PSO没有交叉和变异操作，实现更为简洁；与梯度下降相比，PSO不需要计算导数，能够处理不可导的优化问题。本文将用Java完整实现PSO核心框架，讲解粒子速度与位置更新公式、惯性权重与认知参数的调优策略，并以经典测试函数为案例展示其全局搜索能力。

算法原理

鸟群觅食的数学抽象

在PSO中，每个潜在解被抽象为搜索空间中的一个”粒子”（Particle）。每个粒子具有：

  • 位置向量 xᵢ：表示问题空间中的一个候选解。
  • 速度向量 vᵢ：决定粒子的飞行方向和步长。
  • 个体最优位置 pBestᵢ：该粒子历史上访问过的最优位置。
  • 全局最优位置 gBest：整个群体中所有粒子访问过的最优位置。

速度与位置更新公式

在每次迭代中，粒子的速度和位置按以下公式更新：

vᵢ(t+1) = w·vᵢ(t) + c₁·r₁·(pBestᵢ - xᵢ(t)) + c₂·r₂·(gBest - xᵢ(t))
xᵢ(t+1) = xᵢ(t) + vᵢ(t+1)

其中：
– w：惯性权重（Inertia Weight），控制粒子保持先前速度的程度。
– c₁：认知系数（Cognitive Coefficient），控制粒子向自身历史最优靠拢的程度。
– c₂：社会系数（Social Coefficient），控制粒子向群体最优靠拢的程度。
– r₁, r₂：[0, 1]区间内的随机数，增加搜索的随机性。

参数调优策略

  • 惯性权重 w：
  • 较大值（0.9-1.2）：有利于全局探索，防止陷入局部最优。
  • 较小值（0.4-0.6）：有利于局部开发，加速收敛。

  • 线性递减策略：从0.9线性递减到0.4，兼顾探索与开发。

  • 认知系数 c₁ 与社会系数 c₂：

  • c₁ = c₂ = 2.0 是经典设置。
  • c₁ > c₂：粒子更倾向于个体探索，适合多峰函数。
  • c₁ < c₂：粒子更倾向于群体收敛，适合单峰函数。

核心Java实现

PSO核心框架

import java.util.Arrays;
import java.util.Random;

/**
 * 粒子群优化（PSO）算法实现
 * 支持多维连续函数优化
 */
public class ParticleSwarmOptimization {

    private final int numParticles;    // 粒子数量
    private final int numDimensions;   // 维度数
    private final int maxIterations;   // 最大迭代次数
    private final double w;            // 惯性权重
    private final double c1;           // 认知系数
    private final double c2;           // 社会系数
    private final double[] lowerBound; // 各维度下界
    private final double[] upperBound; // 各维度上界
    private final ObjectiveFunction function; // 目标函数
    private final Random random;

    // 粒子状态
    private double[][] positions;      // 位置
    private double[][] velocities;     // 速度
    private double[][] pBests;         // 个体最优位置
    private double[] pBestValues;      // 个体最优值
    private double[] gBest;            // 全局最优位置
    private double gBestValue;         // 全局最优值

    public ParticleSwarmOptimization(
            int numParticles, int numDimensions, int maxIterations,
            double w, double c1, double c2,
            double[] lowerBound, double[] upperBound,
            ObjectiveFunction function, long seed) {
        this.numParticles = numParticles;
        this.numDimensions = numDimensions;
        this.maxIterations = maxIterations;
        this.w = w;
        this.c1 = c1;
        this.c2 = c2;
        this.lowerBound = lowerBound;
        this.upperBound = upperBound;
        this.function = function;
        this.random = new Random(seed);
    }

    /**
     * 执行PSO优化
     * @return 全局最优位置
     */
    public double[] optimize() {
        initialize();

        for (int iter = 0; iter < maxIterations; iter++) {
            // 线性递减惯性权重
            double currentW = w * (1.0 - (double) iter / maxIterations) + 0.4 * ((double) iter / maxIterations);

            for (int i = 0; i < numParticles; i++) {
                // 评估当前位置
                double fitness = function.evaluate(positions[i]);

                // 更新个体最优
                if (fitness < pBestValues[i]) {
                    pBestValues[i] = fitness;
                    System.arraycopy(positions[i], 0, pBests[i], 0, numDimensions);
                }

                // 更新全局最优
                if (fitness < gBestValue) {
                    gBestValue = fitness;
                    System.arraycopy(positions[i], 0, gBest, 0, numDimensions);
                }
            }

            // 更新速度和位置
            for (int i = 0; i < numParticles; i++) {
                for (int d = 0; d < numDimensions; d++) {
                    double r1 = random.nextDouble();
                    double r2 = random.nextDouble();

                    // 速度更新公式
                    velocities[i][d] = currentW * velocities[i][d]
                            + c1 * r1 * (pBests[i][d] - positions[i][d])
                            + c2 * r2 * (gBest[d] - positions[i][d]);

                    // 限制最大速度（防止飞出搜索空间）
                    double maxVel = (upperBound[d] - lowerBound[d]) * 0.5;
                    velocities[i][d] = Math.max(-maxVel, Math.min(maxVel, velocities[i][d]));

                    // 位置更新
                    positions[i][d] += velocities[i][d];

                    // 边界处理：越界则反弹
                    if (positions[i][d] < lowerBound[d]) {
                        positions[i][d] = lowerBound[d];
                        velocities[i][d] *= -0.5;
                    } else if (positions[i][d] > upperBound[d]) {
                        positions[i][d] = upperBound[d];
                        velocities[i][d] *= -0.5;
                    }
                }
            }

            // 每100轮打印一次进度
            if (iter % 100 == 0 || iter == maxIterations - 1) {
                System.out.printf("迭代 %4d: 全局最优值 = %.10f, 最优位置 = %s%n",
                        iter, gBestValue, arrayToString(gBest));
            }
        }

        return gBest;
    }

    /**
     * 初始化粒子群
     */
    private void initialize() {
        positions = new double[numParticles][numDimensions];
        velocities = new double[numParticles][numDimensions];
        pBests = new double[numParticles][numDimensions];
        pBestValues = new double[numParticles];
        gBest = new double[numDimensions];
        gBestValue = Double.MAX_VALUE;

        for (int i = 0; i < numParticles; i++) {
            for (int d = 0; d < numDimensions; d++) {
                // 随机初始化位置
                positions[i][d] = lowerBound[d] + random.nextDouble() * (upperBound[d] - lowerBound[d]);
                // 随机初始化速度
                velocities[i][d] = (random.nextDouble() - 0.5) * (upperBound[d] - lowerBound[d]) * 0.1;
                pBests[i][d] = positions[i][d];
            }
            pBestValues[i] = Double.MAX_VALUE;
        }
    }

    public double getGBestValue() {
        return gBestValue;
    }

    private String arrayToString(double[] arr) {
        StringBuilder sb = new StringBuilder("[");
        for (int i = 0; i < arr.length; i++) {
            sb.append(String.format("%.6f", arr[i]));
            if (i < arr.length - 1) sb.append(", ");
        }
        sb.append("]");
        return sb.toString();
    }

    /**
     * 目标函数接口
     */
    @FunctionalInterface
    public interface ObjectiveFunction {
        double evaluate(double[] x);
    }
}

经典测试函数

/**
 * 标准优化测试函数库
 * 用于评估PSO等优化算法的性能
 */
public class TestFunctions {

    /**
     * Sphere函数：单峰，全局最小值f(0)=0
     * 用于测试收敛精度
     */
    public static double sphere(double[] x) {
        double sum = 0;
        for (double v : x) sum += v * v;
        return sum;
    }

    /**
     * Rosenbrock函数：单峰但病态，全局最小值f(1,1,...,1)=0
     * 用于测试收敛速度
     */
    public static double rosenbrock(double[] x) {
        double sum = 0;
        for (int i = 0; i < x.length - 1; i++) {
            sum += 100 * Math.pow(x[i + 1] - x[i] * x[i], 2) + Math.pow(1 - x[i], 2);
        }
        return sum;
    }

    /**
     * Rastrigin函数：多峰，全局最小值f(0)=0
     * 用于测试全局搜索能力
     */
    public static double rastrigin(double[] x) {
        double sum = 10 * x.length;
        for (double v : x) {
            sum += v * v - 10 * Math.cos(2 * Math.PI * v);
        }
        return sum;
    }

    /**
     * Ackley函数：多峰，全局最小值f(0)=0
     * 用于测试跳出局部最优的能力
     */
    public static double ackley(double[] x) {
        double sum1 = 0, sum2 = 0;
        for (double v : x) {
            sum1 += v * v;
            sum2 += Math.cos(2 * Math.PI * v);
        }
        return -20 * Math.exp(-0.2 * Math.sqrt(sum1 / x.length))
                - Math.exp(sum2 / x.length) + 20 + Math.E;
    }

    /**
     * Griewank函数：多峰，全局最小值f(0)=0
     */
    public static double griewank(double[] x) {
        double sum = 0, prod = 1;
        for (int i = 0; i < x.length; i++) {
            sum += x[i] * x[i];
            prod *= Math.cos(x[i] / Math.sqrt(i + 1));
        }
        return sum / 4000 - prod + 1;
    }
}

PSO综合测试程序

/**
 * PSO算法综合测试
 * 在多个经典测试函数上验证算法性能
 */
public class PSOTestSuite {

    public static void main(String[] args) {
        int dim = 10;           // 10维问题
        int particles = 50;     // 50个粒子
        int maxIter = 1000;     // 1000次迭代
        double w = 0.9;         // 初始惯性权重
        double c1 = 2.0;        // 认知系数
        double c2 = 2.0;        // 社会系数

        double[] lower = new double[dim];
        double[] upper = new double[dim];

        System.out.println("========== PSO 综合测试 ==========");
        System.out.println("维度: " + dim + ", 粒子数: " + particles + ", 最大迭代: " + maxIter);
        System.out.println("参数: w=" + w + ", c1=" + c1 + ", c2=" + c2);
        System.out.println();

        // 测试1: Sphere函数
        System.out.println("--- 测试1: Sphere函数 ---");
        System.out.println("理论最优值: 0.0 (在x=0处)");
        Arrays.fill(lower, -100);
        Arrays.fill(upper, 100);
        runPSO(particles, dim, maxIter, w, c1, c2, lower, upper, TestFunctions::sphere);

        // 测试2: Rosenbrock函数
        System.out.println("\n--- 测试2: Rosenbrock函数 ---");
        System.out.println("理论最优值: 0.0 (在x=1处)");
        Arrays.fill(lower, -30);
        Arrays.fill(upper, 30);
        runPSO(particles, dim, maxIter, w, c1, c2, lower, upper, TestFunctions::rosenbrock);

        // 测试3: Rastrigin函数
        System.out.println("\n--- 测试3: Rastrigin函数 ---");
        System.out.println("理论最优值: 0.0 (在x=0处)，大量局部最优");
        Arrays.fill(lower, -5.12);
        Arrays.fill(upper, 5.12);
        runPSO(particles, dim, maxIter, w, c1, c2, lower, upper, TestFunctions::rastrigin);

        // 测试4: Ackley函数
        System.out.println("\n--- 测试4: Ackley函数 ---");
        System.out.println("理论最优值: 0.0 (在x=0处)");
        Arrays.fill(lower, -32);
        Arrays.fill(upper, 32);
        runPSO(particles, dim, maxIter, w, c1, c2, lower, upper, TestFunctions::ackley);

        // 测试5: Griewank函数
        System.out.println("\n--- 测试5: Griewank函数 ---");
        System.out.println("理论最优值: 0.0 (在x=0处)");
        Arrays.fill(lower, -600);
        Arrays.fill(upper, 600);
        runPSO(particles, dim, maxIter, w, c1, c2, lower, upper, TestFunctions::griewank);
    }

    private static void runPSO(int particles, int dim, int maxIter, double w, double c1, double c2,
                               double[] lower, double[] upper,
                               ParticleSwarmOptimization.ObjectiveFunction func) {
        // 运行10次取平均
        double bestSum = 0;
        double bestMin = Double.MAX_VALUE;
        double bestMax = 0;

        for (int run = 0; run < 10; run++) {
            ParticleSwarmOptimization pso = new ParticleSwarmOptimization(
                    particles, dim, maxIter, w, c1, c2, lower, upper, func, 42 + run);
            pso.optimize();
            double val = pso.getGBestValue();
            bestSum += val;
            bestMin = Math.min(bestMin, val);
            bestMax = Math.max(bestMax, val);
        }

        System.out.println("10次运行统计:");
        System.out.println("  平均值: " + String.format("%.10f", bestSum / 10));
        System.out.println("  最优值: " + String.format("%.10f", bestMin));
        System.out.println("  最差值: " + String.format("%.10f", bestMax));
    }
}

带约束处理的PSO扩展

/**
 * 带约束处理的PSO实现
 * 使用罚函数法处理不等式约束
 */
public class ConstrainedPSO {

    /**
     * 约束优化问题示例：最小化 f(x,y) = x² + y²
     * 约束: x + y >= 1, x >= 0, y >= 0
     * 理论最优解: (0.5, 0.5), f=0.5
     */
    public static void main(String[] args) {
        int particles = 40;
        int dim = 2;
        int maxIter = 500;

        double[] lower = {0, 0};
        double[] upper = {5, 5};
        double penaltyCoeff = 1000.0; // 罚函数系数

        ParticleSwarmOptimization pso = new ParticleSwarmOptimization(
                particles, dim, maxIter, 0.8, 2.0, 2.0, lower, upper,
                x -> {
                    double f = x[0] * x[0] + x[1] * x[1]; // 目标函数
                    // 约束: x + y >= 1，违反时加惩罚
                    double violation = Math.max(0, 1 - (x[0] + x[1]));
                    return f + penaltyCoeff * violation * violation;
                },
                42
        );

        double[] best = pso.optimize();
        System.out.println("\n带约束优化结果:");
        System.out.println("最优位置: x=" + String.format("%.6f", best[0])
                + ", y=" + String.format("%.6f", best[1]));
        System.out.println("目标函数值: " + String.format("%.6f", best[0] * best[0] + best[1] * best[1]));
        System.out.println("约束满足度: x+y=" + String.format("%.6f", best[0] + best[1]));
    }
}

PSO变体与进阶优化

1. 自适应惯性权重

/**
 * 自适应惯性权重策略
 * 根据群体多样性动态调整w
 */
public static double adaptiveInertiaWeight(double wBase, ParticleSwarmOptimization pso,
                                            double diversityThreshold) {
    double diversity = calculateDiversity(pso);
    if (diversity < diversityThreshold) {
        // 多样性低，增加w以跳出局部最优
        return Math.min(1.2, wBase * 1.1);
    } else {
        // 多样性高，减小w以精细搜索
        return Math.max(0.4, wBase * 0.95);
    }
}

private static double calculateDiversity(ParticleSwarmOptimization pso) {
    // 计算粒子位置的平均欧氏距离
    // 实际实现需要访问粒子位置数据
    return 0.0;
}

2. 拓扑结构变体

拓扑类型 特点 适用场景
全局拓扑 (gBest) 所有粒子共享全局最优 简单问题，快速收敛
环形拓扑 (lBest) 粒子只与邻居共享信息 复杂多峰问题，避免早熟
冯诺依曼拓扑 二维网格邻居结构 平衡探索与开发
随机拓扑 动态随机选择信息来源 高维复杂问题

3. 混合PSO策略

  • PSO + 模拟退火：当群体陷入局部最优时，引入随机扰动。
  • PSO + 遗传算法：定期对粒子进行交叉和变异操作。
  • 多目标PSO (MOPSO)：维护非支配解集，解决多目标优化问题。

复杂度分析

时间复杂度

PSO的单次迭代时间复杂度为 O(N·D)：

  • 评估阶段：N个粒子各评估一次目标函数，通常为O(N·D)。
  • 更新阶段：每个粒子的D维速度和位置更新，为O(N·D)。
  • 总复杂度：O(T·N·D)，其中T是迭代次数，N是粒子数，D是维度数。

相比遗传算法的O(T·N·D·logN)（需要排序选择），PSO的计算效率更高。

空间复杂度

PSO需要存储：
– 粒子位置：O(N·D)
– 粒子速度：O(N·D)
– 个体最优位置：O(N·D)
– 全局最优位置：O(D)

总空间复杂度为 O(N·D)，在内存受限的场景下非常友好。

收敛性分析

PSO的收敛性已被理论证明：

  • 收敛条件：当 w < 1 且 c₁ + c₂ < 4(1+w) 时，粒子期望位置收敛。
  • 早熟收敛：若 c₂ 过大或群体多样性过早丧失，容易陷入局部最优。
  • 参数平衡：w ≈ 0.729, c₁ = c₂ ≈ 1.494 被证明具有良好的收敛性。

PSO的优缺点与适用场景

优点

  • 实现简单：无复杂操作，代码量小，易于理解和部署。
  • 无梯度依赖：仅需目标函数值，适用于不可导、不连续函数。
  • 并行友好：粒子评估完全独立，天然适合GPU并行加速。
  • 全局搜索强：多粒子协同搜索，不易陷入局部最优。
  • 参数少：相比遗传算法，需要调节的参数更少。

缺点

  • 精度有限：对于高精度要求的单峰函数，收敛速度不如梯度方法。
  • 维度诅咒：高维问题（D > 100）中性能下降明显。
  • 参数敏感：w、c₁、c₂的选择对结果影响显著。
  • 离散问题弱：原生PSO针对连续空间，离散问题需要特殊编码。

最佳适用场景

  • 神经网络权重优化：替代反向传播，避免陷入局部最优。
  • 特征选择：在超高维特征空间中寻找最优子集。
  • 调度问题：作业车间调度、任务分配等组合优化。
  • 参数调优：机器学习模型的超参数自动搜索。
  • 路径规划：机器人路径、电路板布线等连续空间搜索。
  • 金融优化：投资组合优化、风险最小化。

总结

粒子群优化是群体智能算法的杰出代表，其简洁的实现和强大的全局搜索能力使其成为解决复杂优化问题的首选工具之一。本文通过Java完整实现了PSO核心框架，包括标准速度-位置更新、线性递减惯性权重、边界处理与约束优化，并在Sphere、Rosenbrock、Rastrigin、Ackley、Griewank五个经典测试函数上验证了算法性能。理解PSO的群体协作机制，不仅能掌握一个实用的优化工具，更能深刻领会”简单规则产生复杂行为”这一涌现计算的核心思想。

Tagged with: