基于C#的粒子群算法（求解双层规划）

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace LBP_ConsoleDemo1
{
    class Program
    {
        static void Main(string[] args)
        {
            
            UInitial();
            DInitial();
            double fitvalue = 0, x = 0, y = 0;
            y = Particle_global_best_position[0];
            for (int j = 0; j < particle_number; j++)
            {
                x = particle_current_position[j, 0];
                fitvalue = UFitness(x, y);
                particle_current_fitness[j] = fitvalue;
                particle_local_best_fitness[j] = particle_current_fitness[j];
            }
            for (int i = 0; i < Generation; i++)
            {
                for (int j = 0; j <=i; j++)
                {
                    DRenewParticle(j);
                }
                URenewParticle(i);
                Console.WriteLine(particle_global_best_fitness + " " + Particle_global_best_fitness);
            }

            Console.WriteLine("测试通过!");
            Console.WriteLine(particle_global_best_fitness + " " + Particle_global_best_fitness);
            Console.ReadKey();
        }


        //***********************************//
        //***********粒子群算法**************//
        //**********************************//


        //粒子群算法容易陷入局部最优

        //指定Particle Swarm Optimizition Algorithm相关参数
        public static double c1 = 2;//"自我认知"系数
        public static double c2 = 2;//"社会认知"系数
        public static double w;//惯性权重
        public static double r = 1;//约束因子
        public static int particle_number = 30;//指定粒子数量
        public static int Generation = 500;//指定迭代次数

        static Random rand = new Random(DateTime.Now.Millisecond);

        //定义粒子属性
        public static int dimention = 1;//指定向量维度(总策略数)

        public static double v_max = 15;//速度最大值
        public static double v_min = 0;//速度最小值
        public static double position_max = 15;//位置最大值
        public static double position_min = 0;//位置最小值

        private static double[,] particle_current_position = new double[particle_number, dimention];//粒子群解集集合
        private static double[,] Particle_current_position = new double[particle_number, dimention];//粒子群解集集合

        private static double[] particle_current_fitness = new double[particle_number];//粒子群各粒子解对应的适应度值
        private static double[] Particle_current_fitness = new double[particle_number];//粒子群各粒子解对应的适应度值

        private static double[,] particle_current_v = new double[particle_number, dimention];//粒子群中各粒子的速度值
        private static double[,] Particle_current_v = new double[particle_number, dimention];//粒子群中各粒子的速度值


        private static double[,] particle_local_best_position = new double[particle_number, dimention];//局部最优粒子群解集集合
        private static double[,] Particle_local_best_position = new double[particle_number, dimention];//局部最优粒子群解集集合


        private static double[] particle_local_best_fitness = new double[particle_number];//局部最优粒子群解对应的局部最优适应度值解集
        private static double[] Particle_local_best_fitness = new double[particle_number];//局部最优粒子群解对应的局部最优适应度值解集


        private static double[] particle_global_best_position = new double[dimention];//粒子群最优解
        private static double[] Particle_global_best_position = new double[dimention];//粒子群最优解


        private static double particle_global_best_fitness;//粒子群全局最优值
        private static double Particle_global_best_fitness;//粒子群全局最优值


        /// <summary>
        /// 适应度值计算
        /// </summary>
        /// <param name="d">解向量</param>
        /// <returns>适应度值</returns>
        private static double UFitness(double x,double y)
        {
            double fitvalue = 0;
            if (-x + y <= 0)
                fitvalue = x * x + (y - 10) * (y - 10);
            else
                fitvalue = x * x + (y - 10) * (y - 10) + 100000 * (-x + y);
            return fitvalue;
        }


        private static double DFitness(double x, double y)
        {
            double fitvalue = 0;
            if (x + y <= 20)
                fitvalue = (x + 2 * y - 30) * (x + 2 * y - 30);
            else
                fitvalue = (x + 2 * y - 30) * (x + 2 * y - 30) + 100000 * (x + y - 20);
            return fitvalue;
        }

        /// <summary>
        /// 粒子群初始化
        /// </summary>
        private static void UInitial()
        {
            //初始化粒子位置、速度
            for (int i = 0; i < particle_number; i++)
            {
                for (int j = 0; j < dimention; j++)
                {
                    particle_current_position[i, j] = rand.Next(15) * rand.NextDouble();
                    particle_local_best_position[i, j] = particle_current_position[i, j];
                    particle_current_v[i, j] = v_min + (v_max - v_min) * rand.NextDouble();
                }
            }
        }

        private static void DInitial()
        {
            //初始化粒子位置、速度
            double x = 0, y = 0;
            for (int i = 0; i < particle_number; i++)
            {
                for (int j = 0; j < dimention; j++)
                {
                    Particle_current_position[i, j] = rand.Next(20) * rand.NextDouble();
                    Particle_local_best_position[i, j] = particle_current_position[i, j];
                    Particle_current_v[i, j] = 20 * rand.NextDouble();
                }
                x = particle_current_position[i, 0];
                y = Particle_current_position[i, 0];
                Particle_current_fitness[i] = DFitness(x, y);
                Particle_local_best_fitness[i] = Particle_current_fitness[i];
            }
            Particle_global_best_fitness = Particle_local_best_fitness[0];
            int index = -1;
            for (int i = 0; i < particle_number; i++)
            {
                if(Particle_local_best_fitness[i]<Particle_global_best_fitness)
                {
                    Particle_global_best_fitness = Particle_local_best_fitness[i];
                    index = i;
                }
            }
            Particle_global_best_position[0] = Particle_local_best_position[index,0];
        }

        private static void URenewParticle(int currentgeneration)
        {
            w = 0.9 - (currentgeneration + 1) * (0.9 - 0.1) / Generation;
            for (int i = 0; i < particle_number; i++)
            {
                double x = 0, y = 0, fitvalue = 0;
                //速度更新
                particle_current_v[i, 0] = w * particle_current_v[i, 0] +
                    c1 * rand.NextDouble() * (particle_local_best_position[i, 0] -
                    particle_current_position[i, 0]) +
                    c2 * rand.NextDouble() * (particle_global_best_position[0] -
                    particle_current_position[i, 0]);
                //位置更新
                particle_current_position[i, 0] += r * particle_current_v[i, 0];
                x = particle_current_position[i, 0];
                y = Particle_global_best_position[0];
                fitvalue = DFitness(x, y);
                if (fitvalue < particle_local_best_fitness[i])
                {
                    particle_local_best_fitness[i] = fitvalue;
                    particle_local_best_position[i, 0] = particle_current_position[i, 0];
                }
            }
            particle_global_best_fitness = particle_local_best_fitness[0];
            for (int i = 0; i < particle_number; i++)
            {
                if (particle_global_best_fitness > particle_local_best_fitness[i])
                {
                    particle_global_best_fitness = particle_local_best_fitness[i];
                    particle_global_best_position[0] = particle_local_best_position[i, 0];
                }
            }
        }

        /// <summary>
        /// 下层种群更新
        /// </summary>
        private static void DRenewParticle(int currentgeneration)
        {
            w = 0.9 - (currentgeneration + 1) * (0.9 - 0.1) / Generation;
            for (int i = 0; i < particle_number; i++)
            {
                double x = 0, y = 0, fitvalue = 0;
                //速度更新
                Particle_current_v[i, 0] = w * Particle_current_v[i, 0] +
                    c1 * rand.NextDouble() * (Particle_local_best_position[i, 0] -
                    Particle_current_position[i, 0]) +
                    c2 * rand.NextDouble() * (Particle_global_best_position[0] -
                    Particle_current_position[i, 0]);
                //位置更新
                Particle_current_position[i, 0] += r * Particle_current_v[i, 0];
                x = particle_current_position[i, 0];
                y = Particle_current_position[i, 0];
                fitvalue = DFitness(x, y);
                if (fitvalue < Particle_local_best_fitness[i])
                {
                    Particle_local_best_fitness[i] = fitvalue;
                    Particle_local_best_position[i, 0] = Particle_current_position[i, 0];
                }
            }
            Particle_global_best_fitness = Particle_local_best_fitness[0];
            for (int i = 0; i < particle_number; i++)
            {
                if (Particle_global_best_fitness > Particle_local_best_fitness[i])
                {
                    Particle_global_best_fitness = Particle_local_best_fitness[i];
                    Particle_global_best_position[0] = Particle_local_best_position[i, 0];
                }
            }
        }
    }
}