一、题目
使用粒子群算法求解函数f(x)的最小值。理论上的最小值是0。
二、原理
粒子群算法利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程,从而获得最优解。
试想一下,如果一群鸟在一片区域中寻找食物,所有的鸟都不知道食物在什么地方,但是每一只鸟都知道自己距离食物有多远,也知道这一群鸟中离食物最近的鸟在什么位置,这样每一只鸟都可以改变当前自己的移动方向,逐渐向离食物最近的鸟所在位置靠近,这样通过不断的搜寻,就能找到食物。
解题思路:
- 假设有100只鸟,初始时这100只鸟分布在不同的位置。这里鸟的位置用30个自变量xi表示。这里的鸟相当于粒子。
- 计算每只鸟相对于食物的距离。这里的距离计算使用设定的目标函数,本题的距离目标函数是f(x)。因为我们就是希望f(x)越小越好,与距离关系对应。
- 找到这100只鸟当中离食物最近的鸟的位置作为此时的全局最优解。
- 根据每一只鸟当前的位置与全局最优解的位置更新当前的速度与位置,向最优解的位置靠近。
- 通过数次迭代,就能找到食物的位置。食物的位置就是f(x)的最小值。
更新公式:
v[i] 表示粒子i的速度。
w 为惯性权重,用于记录当前自身的速度,通常为非负数,调节解空间的搜索范围,为0时则失去自身速度的记忆。
c1、c2 表示学习参数,调节学习的最大步长。
pbest 是当前粒子的最优值。
gbest 为集群中搜索到的最优值,也就是历史最优值。
present[i] 为当前粒子的位置。
三、程序实现
3.1 粒子状态类
import java.util.Random;
public class Particle {
//维数
public int dimension = 30;
//粒子的位置
public double[] X = new double[dimension];
//局部最好位置
public double[] pbest = new double[dimension];
//粒子的速度
public double[] V = new double[dimension];
//适应值
public double fitness;
/**
* 根据当前位置计算适应值
* @return newFitness
*/
public double calculateFitness() {
double newFitness;
double a = 0;
double b = 1;
for (int i = 0; i < dimension; i++) {
a += X[i] * X[i];
b = b * Math.cos(X[i] / Math.pow((i+1), 0.5));
}
newFitness = a / 4000 - b + 1;
return newFitness;
}
/**
* 初始化自己的位置和pbest
*/
public void initialX() {
for(int i = 0;i < dimension; i++) {
X[i] = new Random().nextDouble(); // 随机产生一个0~1的随机小数
X[i] = -600 + 1200 * X[i];
pbest[i] = X[i];
}
}
/**
* 初始化自己的速度
*/
public void initialV() {
for(int i = 0;i < dimension; i++) {
double tmp = new Random().nextDouble(); // 随机产生一个0~1的随机小数
V[i] = tmp * 2 - 4;
}
}
}
3.2 日志类
import java.io.*;
/**
* 日志类,记录结果。
*/
public class Logger {
public static void log(String msg) {
try {
BufferedWriter out = new BufferedWriter(new FileWriter("./PSO/src/log.txt", true));
out.write(msg);
out.flush();
out.close();
} catch (IOException e) {
e.printStackTrace();
}
}
}
3.3 PSO类
import java.util.ArrayList;
import java.util.List;
import java.util.Random;
public class PSO {
private static double[] gbest; // 全局最优位置
private static double gbest_fitness = Double.MAX_VALUE; // 全局最优位置对应的fitness
private static int particle_num = 100; // 粒子数
private static int N = 50000; // 迭代次数
private static int c1 = 2;
private static int c2 = 2;
private static double w = 1.5; // 惯性因子
private static double Vmax = 1; // 最大速度
private static List<Particle> particles = new ArrayList<Particle>();//粒子群
/**
* 主程序入口
*/
public static void main(String[] args) {
initialParticles(); // 初始化
updateGbest(); // 更新gBest
// 循环N次
for (int i = 0; i < N; i++) {
updateV(i); // 更新V
updateX(); // 更新X
updateGbest(); // 更新gBest
// 输出日志
if (i % 20 == 0) {
System.out.println(i + "fitness=" + gbest_fitness);
String msg = i / 20 + " " + gbest_fitness + "\n";
Logger.log(msg);
}
}
}
/**
* 初始化所有粒子
*/
public static void initialParticles() {
for(int i=0; i < particle_num; i++) {
Particle particle = new Particle(); // 创建粒子群对象
particle.initialX(); // 初始化X
particle.initialV(); // 初始化V
particle.fitness = particle.calculateFitness(); // 计算距离
particles.add(particle); // 放到集合中
}
}
/**
* 更新gBest
*/
public static void updateGbest() {
double fitness = Double.MAX_VALUE;
int index = 0;
for(int i = 0; i < particle_num; i++) { // 找到群体中适应值最小的粒子
if(particles.get(i).fitness<fitness) {
index = i;
fitness = particles.get(i).fitness;
}
}
if(fitness < gbest_fitness) { // 如果个体适应值小于全局适应值,更新全局的最优值为个体最优值
gbest = particles.get(index).pbest.clone();
gbest_fitness = fitness;
}
}
/**
* 更新每个粒子的速度
*/
public static void updateV(int n) {
if (n % 10000 == 0) {
Vmax = Vmax * 0.1;
}
for(Particle particle : particles) {
for(int i = 0;i < particle.dimension; i++) {
double v = w * particle.V[i]+c1*rand()*(particle.pbest[i]-particle.X[i])+c2*rand()*(gbest[i]-particle.X[i]);
if(v > Vmax) // 判断速度是否超过最大的速度
v = Vmax;
else if(v < - Vmax) // 比最大速度的相反数小
v = - Vmax;
particle.V[i] = v;//更新Vi
}
}
}
/**
* 更新每个粒子的位置和pbest
*/
public static void updateX() {
for(Particle particle:particles) {
for(int i=0;i<particle.dimension;i++) {
particle.X[i] = particle.X[i] + particle.V[i];
}
double newFitness = particle.calculateFitness();//新的适应值
//如果新的适应值比原来的小则跟新fitness和pbest
if(newFitness < particle.fitness) {
particle.pbest = particle.X.clone();
particle.fitness = newFitness;
}
}
}
/**
* 返回一个0~1的随机数
* @return
*/
public static double rand() {
return new Random().nextDouble();
}
}
四、实验结果
4.1 Python可视化
这里为了显示的美观性,只选取了其中一部分数据进行可视化。
import os
import matplotlib.pyplot as plt
import matplotlib; matplotlib.use('TkAgg')
# 读取log.txt文件
with open(os.path.join("log.txt"), 'r') as fd:
dataList = fd.readlines()
dataList = [l[:-1] for l in dataList] # 把/n去掉
# 将数据放入list中
data = []
iters = []
for i in range(len(dataList[900:950])):
sample = dataList[i].split()
data.append(float(sample[1]))
iters.append(int(sample[0]))
# 绘图
plt.title("Result GA")
plt.plot(iters, data, color='blue')
plt.xlabel("iters")
plt.ylabel("result")
plt.savefig("Result_GA")
plt.show()
4.2 实验结果
横坐标是迭代次数,纵坐标是距离结果。
迭代50000次的结果。理论上继续迭代可以使值更接近于0。
五、总结
本次实验在调参的时候遇到了问题,无论迭代次数有多少次,目标距离很难有进一步的缩小,即无法跳出局部最优解,经过大量的调试,得到了其中一种解决方案。
方法是改变设定的Vmax值,该值是用来限制V的大小,V是用来更新参数X的,X是用来求距离的。一共50000次迭代,每10000次,Vmax缩小为原来的十分之一。
通过该方法,解决了距离无法缩小的问题。即无法跳出局部最优解的问题。