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++){//计算mbestfor (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;elseswarm[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()中的结果输出，运行过程及结果如下图所示。