【算法】通过粒子群优化算法pso采用python和c#进行多变量函数的最小值求解

---

需求

求解f=sin（x1）+…sin(x50）最小值，其中x取值范围各不同，上下限为【0,3】之间任意小数。通过粒子群优化算法pos方法，采用python和c#进行求解

粒子群优化算法PSO

粒子群优化算法（Particle Swarm Optimization，简称PSO）是一种启发式优化算法，类似于群体智能算法中的“鸟群”行为，通过模拟鸟群中个体之间的信息交流和协作，来寻找最优解。

基本思想

将待求解问题的解看作粒子，在解空间中随机分布，每个粒子根据自身的历史经验和整个群体的最优解进行搜索，并调整自身的位置和速度，以找到更优的解。整个算法的运行过程可以分为初始化、迭代更新和终止条件三个阶段。

• 初始化阶段，需要设置粒子群的初始位置和速度，并随机生成每个粒子的个体和全局最优解。
• 迭代更新阶段，每个粒子根据自身的历史最优解和整个群体的最优解，通过更新速度和位置，来寻找更优解。具体步骤如下：

1. 计算每个粒子的适应度值，并更新个体最优解；
2. 更新全局最优解；
3. 根据当前位置和速度更新粒子的速度和位置；
4. 判断是否满足终止条件，如果满足则算法结束，否则进入下一次迭代。

• 终止条件可以根据问题的实际情况设置，例如达到最大迭代次数、适应度值收敛等。

优点

1. 算法简单
   PSO算法的原理和实现相对简单，易于理解和实施。
2. 全局搜索能力
   PSO算法具有较好的全局搜索能力，在多维连续空间中能够快速找到全局最优解。
3. 并行计算
   PSO算法是一种群体智能算法，可以进行并行计算，适合于分布式计算和并行处理。
4. 鲁棒性
   PSO算法对问题的初始点位置不敏感，能够在大部分情况下得到较好的结果。

缺点

1. 参数选择困难
   PSO算法涉及到较多的参数选择，如粒子数量、惯性权重、加速因子等，参数选择不合适会影响算法的性能。
2. 陷入局部最优
   PSO算法容易陷入局部最优解，特别是在复杂的多峰函数中，全局最优解的搜索能力有限。
3. 不适合离散问题
   PSO算法主要适用于连续优化问题，对于离散问题的处理较为困难。
4. 收敛速度较慢
   PSO算法的收敛速度较慢，特别是在高维问题中，需要较长的迭代次数才能达到较好的结果。

综上所述，粒子群优化算法具有简单、全局搜索能力强等优点，但在参数选择、陷入局部最优等方面存在一定的缺点。在实际应用中，需要根据具体问题的特点和要求来选择合适的优化算法。

python方法

采用PyCharm Community Edition 2021.2.3编译器进行编写，通过pso方法进行函数求解

import csv
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
from datetime import datetime
from sklearn.metrics import explained_variance_score
from sklearn.svm import SVR
from sklearn.metrics import mean_squared_error
from sklearn.metrics import accuracy_score
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.metrics import explained_variance_score
from sklearn import metrics
from sklearn.metrics import mean_absolute_error # 平方绝对误差
import random
import pandas as pd
import math


class PSO:
    def __init__(self, parameters):
        """
        particle swarm optimization
        parameter: a list type, like [NGEN, pop_size, var_num_min, var_num_max]
        """
        # 初始化
        self.NGEN = parameters[0]  # 迭代的代数
        self.pop_size = parameters[1]  # 种群大小
        self.var_num = len(parameters[2])  # 变量个数
        self.bound = []  # 变量的约束范围
        self.bound.append(parameters[2])
        self.bound.append(parameters[3])

        self.pop_x = np.zeros((self.pop_size, self.var_num))  # 所有粒子的位置
        self.pop_v = np.zeros((self.pop_size, self.var_num))  # 所有粒子的速度
        self.p_best = np.zeros((self.pop_size, self.var_num))  # 每个粒子最优的位置
        self.g_best = np.zeros((1, self.var_num))  # 全局最优的位置

        # 初始化第0代初始全局最优解
        # temp = -1
        temp = 50
        for i in range(self.pop_size):
            for j in range(self.var_num):
                self.pop_x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j])
                self.pop_v[i][j] = random.uniform(0, 1)
            self.p_best[i] = self.pop_x[i]  # 储存最优的个体
            fit = self.fitness(self.p_best[i])
            if fit < temp:
                self.g_best = self.p_best[i]
                temp = fit

    def fitness(self, ind_var):

        a = 0
        for i in range(0, 50):
            a = a + np.sin(ind_var[i])
        return a

    def update_operator(self, pop_size):
        """
        更新算子：更新下一时刻的位置和速度
        """
        c1 = 2  # 学习因子，一般为2
        c2 = 2
        w = 0.4  # 自身权重因子
        for i in range(pop_size):
            # 更新速度
            self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * (
                    self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i])
            # 更新位置
            self.pop_x[i] = self.pop_x[i] + self.pop_v[i]
            # 越界保护
            for j in range(self.var_num):
                if self.pop_x[i][j] < self.bound[0][j]:
                    self.pop_x[i][j] = self.bound[0][j]
                if self.pop_x[i][j] > self.bound[1][j]:
                    self.pop_x[i][j] = self.bound[1][j]
            # 更新p_best和g_best
            if self.fitness(self.pop_x[i]) < self.fitness(self.p_best[i]):
                self.p_best[i] = self.pop_x[i]
            if self.fitness(self.pop_x[i]) < self.fitness(self.g_best):
                self.g_best = self.pop_x[i]

    def main(self):
        popobj = []
        # self.ng_best = np.zeros((1, self.var_num))[0]
        for gen in range(self.NGEN):
            self.update_operator(self.pop_size)
            # popobj.append(self.fitness(self.g_best))
            # print('############ Generation {} ############'.format(str(gen + 1)))
            # if self.fitness(self.g_best) < self.fitness(self.ng_best):
            #     self.ng_best = self.g_best.copy()
        print('最好的位置：{}'.format(self.g_best))
        print('最小的函数值：{}'.format(self.fitness(self.g_best)))
        print("---- End of (successful) Searching ----")


if __name__ == '__main__':
    NGEN = 200
    popsize = 20
    low = []
    up = []
    for i in range(0,50):
        a=round(random.uniform(0,3),2)
        b=round(random.uniform(0,3),2)
        if(a >b):
            low.append(b)
            up.append(a)
        else:
            low.append(a)
            up.append(b)
    print(low)
    print(up)
    parameters = [NGEN, popsize, low, up]
    pso = PSO(parameters)
    pso.main()

c#方法

采用Visual Studio Code编译器进行编写，通过pos方法进行函数求解

// See https://aka.ms/new-console-template for more information
// Console.WriteLine("Hello, Worlddddd!");
using System;

namespace test
{
    class Program
    {
        /* 已知函数 𝑦=𝑓(𝑥1,𝑥2)=𝑥1²+𝑥2² ，其中−10≤𝑥1,𝑥2≤10
        请用粒子群优化算法求解y的最小值 */

        static int N = 100;//初始化种群大小
        static int x_number = 50;//变量个数
        int minValue = 0;
        int maxValue = 3;//函数自变量取值范围
        double[] min_xreal = new double[x_number];
        double[] max_xreal = new double[x_number];
        double v_max = 0.5;
        double v_min = -0.5;//速度范围
        double w = 0.5;//惯量权重
        double c1 = 2;
        double c2 = 2;//加速系数
        Random r = new Random();//随机数
        double[,] particle = new double[N, x_number];//初始化3个粒子
        double[,] current_v = new double[N, x_number];//粒子当前速度
        double[,] pBest = new double[N, x_number];//局部最优解
        double[] gBest = new double[x_number];//全局最优解
        double[] current_Fitness = new double[N];//粒子当前适应度（函数值）
        double[] particle_local_best_Fitness = new double[N];//粒子局部最优适应度（函数值）
        double particle_global_Fitness;//粒子全局最优适应度（函数值）
        public void Initial()//粒子初始化
        {
            int i, j;

            //确定变量范围
            for (int a = 0; a < x_number; a++)
            {
                double r_temp = r.NextDouble() * (maxValue - minValue) + minValue;
                double r_temp1 = r.NextDouble() * (maxValue - minValue) + minValue;
                if (r_temp > r_temp1)
                {
                    min_xreal[a] = r_temp1;
                    max_xreal[a] = r_temp;
                }
                else
                {
                    min_xreal[a] = r_temp;
                    max_xreal[a] = r_temp1;
                }
            }

            for (i = 0; i < N; i++)
            {
                for (j = 0; j < x_number; j++)
                {
                    double min, max = 0;
                    min = min_xreal[j];
                    max = max_xreal[j];
                    particle[i, j] = min + (max - min) * r.NextDouble();//初始化群体
                    pBest[i, j] = particle[i, j];//将当前最优结果写入局部最优
                    current_v[i, j] = v_min + 2 * v_max * r.NextDouble(); //初始化粒子的速度
                }

            }

            for (i = 0; i < N; i++) //计算每个粒子的适应度
            {
                double[] tmp = new double[x_number];
                for (int h = 0; h < x_number; h++)
                {
                    tmp[h] = particle[i, h];
                }
                current_Fitness[i] = Fitness(tmp);
                particle_local_best_Fitness[i] = current_Fitness[i];
            }

            particle_global_Fitness = particle_local_best_Fitness[0];//找出全局最优的适应度(函数值)
            j = 0;
            for (i = 0; i < N; i++)
            {
                if (particle_local_best_Fitness[i] < particle_global_Fitness)
                {
                    particle_global_Fitness = particle_local_best_Fitness[i];
                    j = i;
                }
            }

            for (i = 0; i < x_number; i++) //更新全局最优向量
            {
                gBest[i] = pBest[j, i];
            }
            Console.WriteLine("初始化粒子完成!");
            Console.WriteLine("正在进行迭代...");
        }

        public void Renew_location()//迭代更新粒子所在的位置和速度
        {
            for (int i = 0; i < N; i++)
            {
                for (int j = 0; j < x_number; j++)
                {
                    double r1 = r.NextDouble();
                    double r2 = r.NextDouble();
                    current_v[i, j] += (w * current_v[i, j] + c1 * r1 * (pBest[i, j] - particle[i, j]) +
                        c2 * r2 * (gBest[j] - particle[i, j]));//更新速度
                    particle[i, j] += current_v[i, j];//更新位置
                    if (particle[i, j] > max_xreal[j])//越界的位置，合法性调整(定义域范围内)
                    {
                        particle[i, j] = max_xreal[j];
                    }
                    if (particle[i, j] < min_xreal[j])
                    {
                        particle[i, j] = min_xreal[j];
                    }
                }
            }

        }
        static double Fitness(double[] x1)            //返回y=f(x1,x2)=x1²+x2²的函数值
        {
            double result = 0;
            for (int i = 0; i < x_number; i++)
            {
                result = result + Math.Sin(x1[i]);
            }
            return result;
        }


        public void Renew_Fitness()//评估例子的适应度函数值
        {
            int j = -1;
            for (int i = 0; i < N; i++)
            {
                double[] tmp = new double[x_number];
                for (int h = 0; h < x_number; h++)
                {
                    tmp[h] = particle[i, h];
                }
                if (Fitness(tmp) < current_Fitness[i])//更新局部最优解pBest
                {
                    for (int a = 0; a < x_number; a++) //更新局部最优向量
                    {
                        pBest[i, a] = particle[i, a];
                    }

                }
                double[] tmp1 = new double[x_number];
                for (int h1 = 0; h1 < x_number; h1++)
                {
                    tmp1[h1] = pBest[i, h1];
                }
                particle_local_best_Fitness[i] = Fitness(tmp1);//更新局部最优适应值

                if (particle_local_best_Fitness[i] < particle_global_Fitness)//更新全局最优适应度
                {
                    particle_global_Fitness = particle_local_best_Fitness[i];
                    j = i;
                    for (int a = 0; a < x_number; a++) //更新全局最优向量
                    {
                        gBest[a] = pBest[j, a];
                    }
                }

            }


        }
        public static void PrintArray(double[] a)
        {
            foreach (var x in a)
            {
                Console.Write("{0}", x);
                Console.Write("/");
            }
        }
        static void Main(string[] args)
        {
            Program PSO = new Program();
            PSO.Initial();
            for (int i = 0; i < 100; i++)//迭代100次
            {
                PSO.Renew_location();
                PSO.Renew_Fitness();
            }
            Console.WriteLine("下限取值");
            PrintArray(PSO.min_xreal);
            Console.WriteLine("/");
            Console.WriteLine("上限取值");
            PrintArray(PSO.max_xreal);
            Console.WriteLine("/");
            Console.WriteLine("取值");
            PrintArray(PSO.gBest);
            Console.WriteLine("迭代完成，最小值y={0}", PSO.particle_global_Fitness);
            Console.ReadLine();
        }
    }
}