PSO算法

PSO算法

作为一种演化算法, PSO(Particle Swarm Optimization)因为它的易实现和参数少的特性被广泛使用.

虽然我不能理解怎么能称得上是人工智能, 但是似乎沾一点边, 和SOM有一点点像.

逻辑

个人与群体

使用粒子群优化, 很显然, 需要操控粒子在每次迭代中不断接近最优值.

我们一开始随机生成选择的粒子, 需要有不断向全局最优进发的趋势.

所以, 在每次迭代时PSO会将粒子位置进行更新, 这就需要计算下次更新应前进的位置.

首先, 对每个粒子来说, 它总倾向于向它所经历过的个人最优点(pbest)进发, 但是后果就是很容易陷入局部最优陷阱.

那么我们就需要对粒子的前进进行更多的调控, 这里我们引入当前的全局最优点(gbest), 并让所有粒子都拥有向当前全局最优进发的趋势. 同时, 对个人最优和全局最优加权.

到此, 迭代的位移公式可以写成

P_{coord}^{(i)} = P_{coord}^{(i)} + w_1(P_{pbest}^{(i)}-P_{coord}^{(i)}) + w_2(P_{gbest}^{(i)}-p_{coord}^{(i)})

有了个人最优解和全局最优解, 我们就可以寻找空间中的最优解了.

惯性

我们期望的是, 在经历有限次迭代后, 粒子能收敛于全局最优点, 经过实测, 在极其有限的迭代次数中, 光靠个人和群体位移是无法让所有点都收敛于全局最优的, 所以我们引入一个惯性权重, 将粒子前一时刻的位移加权后加入下次的位移(称之为惯性十分贴切), 这样能够有效的解决粒子在两个或多个点之间反复横跳而无法收敛的问题.

到此, 迭代的位移公式则为:

P_{coord}^{(i)} = w_0P_{coord}^{(i)} + w_1(P_{pbest}^{(i)}-P_{coord}^{(i)}) + w_2(P_{gbest}^{(i)}-p_{coord}^{(i)})

更新时到底在做什么?

首先, 在更新位置时, 我们需要对边界进行限制, 可能的情况下, 速度也要进行适当限制(不确定效果), 所以在更新前要对公式进行限幅.

P_{coord}^{(i)} =\LARGE{Max\{}\normalsize{Max\{w_0P_{coord}^{(i)} + w_1(P_{pbest}^{(i)}-P_{coord}^{(i)}) + w_2(P_{gbest}^{(i)}-p_{coord}^{(i)}),velocitylimits\}\LARGE{,boundaries\}}}

在代码实现中, 使用循环和辅助变量即可完成此操作.

其次, 合理地使用辅助变量和函数进行额外操作也十分重要, 寻找全局最大值时可以自己写简单的遍历, 也可以使用数据结构中自带的方法求最大值.

真代码

二维下的未进行权重优化的实现:

using System;
using System.Collections.Generic;
using System.Reflection;
using System.Reflection.Emit;
using System.Text;

namespace PSO_CS
{
    public class PSO
    {
        struct Particle {   //为粒子提供数据结构
            public double x, y,bestx,besty;
            public double velocityX,velocityY;
            public double fitness;
            public double bestfitness;
        }

        int numberOfParticle;   //粒子数
        double c1,c2,w; //自我学习因子, 群体学习因子, 惯性权重
        double minX, minY, maxX, maxY;  //边界
        int max_iteration;  //最大循环次数
        double velocityLimitX,velocityLimitY;   //速度限制
        private double globalBestFitness = Double.MinValue ;   //全局最大fitness
        private double globalBestX, globalBestY;    //全局最大位置
        private Particle[] particles;   //粒子实例
        public int NumberOfParticle { get => numberOfParticle; }
        public double C1 { get => c1; }
        public double C2 { get => c2; }
        public double W { get => w; }
        public double MinX { get => minX; }
        public double MinY { get => minY; }
        public double MaxX { get => maxX; }
        public double MaxY { get => maxY; }
        public int Max_iteration { get => max_iteration; }
        public double VelocityLimitX { get => velocityLimitX; }
        public double VelocityLimitY { get => velocityLimitY; }

        public PSO()    //构造器
        {
            InitParicle(); 
        }

        private double FitnessFunction(double x,double y) { //对象
                return Math.Cos(x)*Math.Sin(x*y)-Math.Sin(y)*Math.Cos(x*y)+x*Math.Cos(y+Math.PI*x)+x*Math.Cos(y)/5;
        }

        private void UpdatedFitnessFunction() { //更新粒子fitness
            for (int i = 0; i < numberOfParticle; i++) {
                particles[i].fitness = FitnessFunction(particles[i].x, particles[i].y);
            }
        }

        private void RandomizedParticle() { //随机粒子位置
            Random random = new Random();
            for (int i = 0; i < numberOfParticle; i++) {
                particles[i].x = (maxX - minX) * random.NextDouble() + minX;
                particles[i].y = (maxY - minX) * random.NextDouble() + minY;
                particles[i].velocityX = particles[i].velocityY = 0;
            }
        }

        private void InitParicle() {    //初始化
            Console.Write("请输入粒子数:");
            numberOfParticle = int.Parse(Console.ReadLine());
            Console.Write("请输入x下界、上界，y下界、上界(分为四次):");
            minX = double.Parse(Console.ReadLine());
            maxX = double.Parse(Console.ReadLine());
            minY = double.Parse(Console.ReadLine());
            maxY = double.Parse(Console.ReadLine());
            velocityLimitX = (maxX - minX) / 2; //草率的速度限制
            velocityLimitY = (maxX - minX) / 2;
            Console.Write("请输入迭代次数:");
            max_iteration = int.Parse(Console.ReadLine());
            w = 0.8;
            c1 = 0.5;
            c2 = 0.5;
            particles = new Particle[numberOfParticle]; //实例初始化
            RandomizedParticle();   //随机粒子位置
            for (int i = 0; i < numberOfParticle; i++) {    //初始化参数
                particles[i].bestfitness = 0;
                particles[i].bestx = 0;
                particles[i].besty = 0;
            }
            UpdatedFitnessFunction();   //更新权重(第0次)
        }
        /// <summary>
        /// 寻找全局最优位置和值
        /// </summary>
        /// <param name="fitness"></param>
        /// <returns>该方法返回x,y,fitness元组</returns>
        private Tuple<double,double,double> MaxFitness(ref double[,] fitness) {
            double max = Double.MinValue;
            double x=double.NaN,y= double.NaN;
            for (int i = 0; i < fitness.Length/3; i++) {
                if (fitness[i, 2] > max) {
                    max = fitness[i, 2];
                    x = fitness[i, 0];
                    y = fitness[i, 1];
                }
            }
            return new Tuple<double,double,double>(x,y,max);
        }

        public void Run()   //算法实现
        {
            int iteration = 0;
            Random random = new Random();
            while (iteration < max_iteration)   //小于迭代数
            {
                //x,y,fitness的辅助数组
                double[,] fitness = new double[numberOfParticle,3]; 

                for (int i = 0; i < numberOfParticle; i++)  //复制权重和位置
                {
                    fitness[i, 0] = particles[i].x;
                    fitness[i, 1] = particles[i].y;
                    fitness[i,2] = FitnessFunction(particles[i].x, particles[i].y);
                }
                UpdatedFitnessFunction();   //更新粒子权重
                for (int i = 0; i < numberOfParticle; i++) {//更新个体最优值
                    if (particles[i].bestfitness < fitness[i,2]){
                        particles[i].bestx = fitness[i, 0];
                        particles[i].besty = fitness[i, 1];
                        particles[i].bestfitness = fitness[i,2];
                    }
                }

                double max, x, y;
                max = MaxFitness(ref fitness).Item3;
                x = MaxFitness(ref fitness).Item1;
                y= MaxFitness(ref fitness).Item2;
                //找最大fitness
                if (globalBestFitness<max) {
                    globalBestFitness = max;
                    globalBestX = x;
                    globalBestY = y;
                }

                for (int i = 0; i < numberOfParticle; i++)  //更新位置
                {
                    double vx = particles[i].velocityX * w +
                        c1 * random.NextDouble() * (particles[i].bestx - particles[i].x) +
                        c2 * random.NextDouble() * (globalBestX - particles[i].x);
                    double vy = particles[i].velocityY * w +
                        c1 * random.NextDouble() * (particles[i].bestx - particles[i].x) +
                        c2 * random.NextDouble() * (globalBestY - particles[i].x);
                    particles[i].velocityX = vx < velocityLimitX ? vx > -velocityLimitX ? vx : -velocityLimitX : velocityLimitX;
                    particles[i].velocityY = vy < velocityLimitY ? vy > -velocityLimitY ? vy : -velocityLimitY : velocityLimitY;
                    double tempX = particles[i].velocityX + particles[i].x;
                    double tempY = particles[i].velocityY + particles[i].y;
                    particles[i].x = tempX < maxX ? tempX > minX ? tempX : minX :maxX;
                    particles[i].y = tempY < maxY ? tempY > minX ?  tempY : minY:maxY;
                }
                iteration++;
            }
            Console.WriteLine("迭代:" + max_iteration);
            Console.WriteLine("结果: (" + globalBestX + ',' + globalBestY + ") with value " + globalBestFitness);
        }
    }
}

实例化对象再调用Run(),

输入参数后

差距很小但仍有待优化.

缺陷

其实这个算法的最大缺陷在于不精确性. 不过我还在研究中, 请等待下一次更新吧.

优化

关于惯性权重, 线性变化惯性权重据说能让算法更快和更准确的收敛.

还有模糊权重、粒子淘汰等等, 其实框架有了还挺简单.