QPSO Algorithm
C#语言.NetFramwork4.6.1平台实现(需了解QPSO算法原理,可参考清华大学孙俊教授编写的教材《量子行为粒子群优化原理及其应用》)
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
namespace Quantum_Behavioral_Particle_Swarm_Algorithm
{
public struct Particle //定义粒子的结构体
{
public double[] x;//粒子当前位置
public double[] pbest;//粒子的个体最好位置
public double bestfitness;//粒子个体最好位置适应值
public double fitness;//粒子的当前最好位置适应值
}
public class QPSO_Test
{
//全局变量定义
readonly int popsize = 20;//群体规模
readonly static int dimension = 30;//问题维数
readonly static int runno = 50;//运行次数
readonly int Maxiter = 2000;//最大迭代次数
double irange_l = 15;//位置初始化下界
double irange_r = 30;//位置初始化上界
double xmin = -100;//搜索上界
double xmax = 100;//搜索下界
double[] gbest = new double[dimension];//全局最优解
double gbestfitness=0;//全局最优值
Random rand = new Random();//随机数产生器
public static double[] data = new double[runno];//存储每次运行后的结果
public static double[] gdata = new double[runno];//种群迭代结果存储
/// <summary>
/// 标准测试函数Sphere函数定义
/// </summary>
/// <param name="x">可行解</param>
/// <returns>函数值</returns>
double Fsphere(double[] x)
{
int i;
double value;
value = 0;
for (i = 0; i < dimension; i++)
{
value += Math.Pow(x[i], 2);
}
return value;
}
/// <summary>
/// 标准测试函数Rosenbrock函数的定义
/// </summary>
/// <param name="x">可行解</param>
/// <returns>函数值</returns>
double Frosenbrock(double[] x)
{
int i;
double value;
value = 0;
for (i = 0; i < dimension-1; i++)
{
value += 100 * Math.Pow((x[i + 1] - x[i] * x[i]), 2) + Math.Pow((x[i] - 1), 2);
}
return value;
}
/// <summary>
/// 程序初始化定义
/// </summary>
/// <param name="population">粒子群种群</param>
void Initiate(Particle[] population)
{
int i, j;
for (i = 0; i < population.Length; i++)
{
population[i].x = new double[dimension];
population[i].pbest = new double[dimension];
for (j = 0; j < dimension; j++)
{
//粒子当前位置的初始化
population[i].x[j] = rand.NextDouble() * (irange_r - irange_l) + irange_l;
//粒子当前最好位置的初始化
population[i].pbest[j] = population[i].x[j];
}
//粒子当前位置的适应值
population[i].fitness = Frosenbrock(population[i].x);
//粒子个体最好位置的适应值
population[i].bestfitness = population[i].fitness;
}
}
/// <summary>
/// 更新粒子群全局最好位置
/// </summary>
/// <param name="population">粒子群</param>
/// <returns>最好解索引下标</returns>
int Globalbest(Particle[] population)
{
int i, flag;
double s = 0;
s = population[0].fitness;
flag = 0;
for (i = 1; i < population.Length; i++)
{
if (population[i].fitness < s)
{
s = population[i].fitness;
flag = i;
}
}
return flag;
}
/// <summary>
/// 主程序
/// </summary>
public void RunAlgorithm()
{
Particle[] swarm = new Particle[popsize];//定义粒子群
//其他变量
int i, k, t, g, run;//run:运行次数循环下标 k:维数循环下标 t:迭代次数循环下标 i:种群循环下标 g:最好位置粒子下标
double a, tmp, fi1, fi2, u, v, z, b, p;//a:收缩-扩张系数值 tmp:临时局部变量 p:吸引势量子 fi1、fi2、u、z:随机数 b:特征长度 v:随机对数
double[] mbest = new double[dimension];//每维平均位置
double[] gbest = new double[dimension];//全局最好粒子
//循环迭代寻优
for (run = 0; run < runno; run++)
{
Initiate(swarm);//初始化群体
//求全局最好位置粒子下标
g = Globalbest(swarm);
for (k = 0; k < dimension; k++)
{
gbest[k] = swarm[g].pbest[k];
}
gbestfitness = swarm[g].bestfitness;//初始化最优值的变量
//算法迭代开始
int j = 0;
for (t = 0; t < Maxiter; t++)
{
//计算mbest
for (k = 0; k < dimension; k++)
{
tmp = 0;
for (i = 0; i < popsize; i++)
{
tmp += swarm[i].pbest[k];
}
mbest[k] = 1.0 * tmp / popsize;
}
//收缩-扩张系数值的计算
a = (1.0 - 0.5) * (Maxiter - t) / Maxiter + 0.5;
//群体更新
for (i = 0; i < popsize; i++)
{
for (k = 0; k < dimension; k++)
{
fi1 = rand.NextDouble();
fi2 = rand.NextDouble();
p = 1.0 * (fi1 * swarm[i].pbest[k] + fi2 * gbest[k]) / (fi1 + fi2);
u = rand.NextDouble();
b = a * Math.Abs(mbest[k] - swarm[i].x[k]);
v = Math.Log(1.0 / u);
z = rand.NextDouble();
//更新粒子位置
if (z < 0.5)
swarm[i].x[k] = p + b * v;
else
swarm[i].x[k] = p - b * v;
//粒子位置限制在搜索范围内
if (swarm[i].x[k] < xmin)
swarm[i].x[k] = xmin;
if (swarm[i].x[k] > xmax)
swarm[i].x[k] = xmax;
}
swarm[i].fitness = Frosenbrock(swarm[i].x);//粒子当前位置适应值
if (swarm[i].fitness < swarm[i].bestfitness)//更新粒子个体最好位置
{
for (k = 0; k < dimension; k++)
{
swarm[i].pbest[k] = swarm[i].x[k];
}
swarm[i].bestfitness = swarm[i].fitness;
}
if (swarm[i].bestfitness < gbestfitness)//更新粒子全局最好位置
{
for (k = 0; k < dimension; k++)
{
gbest[k] = swarm[i].pbest[k];
}
gbestfitness = swarm[i].bestfitness;
}
}
if ((t+1) % 40 == 0)
{
gdata[j++] = gbestfitness;
}
}
data[run] = gbestfitness;
}
}
/// <summary>
/// 算法运行结果
/// </summary>
/// <returns>最优值</returns>
public double[] Value()
{
return data;
}
/// <summary>
/// 种群迭代结果
/// </summary>
/// <returns>迭代值</returns>
public double[] GValue()
{
return gdata;
}
}
}
运行结果图
运用Chart控件可将方法Value()与方法GValue()中的结果输出,运行过程及结果如下图所示。