QPSO Algorithm

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()中的结果输出，运行过程及结果如下图所示。