高斯协同粒子群优化算法的实现

实验报告题目：多尺度协同变异的粒子群优化算法的实现学院：信息科学与技术学院专业：计算机技术年级：2016级研究生课程：计算智能_作业4姓名：学号：备注：一、实验内容1.问题描述粒子群算法(Particle Swarm Optimizer，PSO)是由Kennedy 和Eberhart 博士提出的一种基于群体智能的优化算法,其基本思想是受到他们早期对许多鸟类的群体行为进行建模与仿真研究的启发.粒子群算法的优势在于其简单容易实现,没有很多参数需要调整,是非线性连续优化问题、组合优化问题和混合整数非线性优化问题的有效优化工具.由于PSO 算法概念简单、实现容易，在函数优化和神经网络权值训练等方面都有很好的表现.由于其又具有深刻的智能背景，既适合科学研究，又特别适于工程应用,因此在近年来得到了飞速的发展.其应用涉及系统控制、人工智能、模式识别、计算机、通信工程等各个领域.粒子群优化算法问世以来受到了广泛的重视,经过十几年的研究,其不论在应用方面还是在优化性能方面都得到了很大的发展.然而,研究人员仍然无法解决一直困扰该算法的两个难题:早熟和收敛速度慢.在避免早熟收敛的研究中,一些学者提出了控制种群多样性来提高性能的办法.其思想主要分为两个方面:（1）从PSO 参数以及其拓扑结构研究的角度.（2）从与遗传算法相结合的角度.对于第一个角度,一些学者提出解决微粒间的冲突和聚集,分析微粒和种群最优位置的距离,种群随机多代初始化的思想.而对于第二个角度,主要是通过引入遗传算法中的变异操作实现的,其中文献通过对微粒引入一个小概率随机变异操作来增加种群多样性（dissipative PSO 简称DPSO）,这样做虽然可以更好的搜索粒子群的解空间,但如变异率控制不当则极容易导致原始种群的混乱,致使在算法的后期达不到精确的局部搜索目的;在此基础上,进一步提高搜索性能,文献给出了一种令变异操作随时间变化的自适应层次PSO 算法(self-organizinghierarchical PSO，简称HPSO),为了使其能够适应变异操作,还确定了自适应参数的选择方式.但由于该算法没有考虑速度公式中惯性因素的影响,因此将发生变异的条件设置为微粒的速度为0,从而使得该算法仍然未能迅速有效的逃离局部极小值.文献为了提高种群多样性和搜索速度间的平衡,对阈值设定进行了改进并结合群体自动分家机制,提出一种自适应逃逸的微粒群算法(Self-Adaptive Escape PSO,简称AEPSO).然而上述算法在实现变异的时候,均采用单一的均匀变异机制,其逃逸能力大大减弱,进而使改进的PSO 算法在有限的迭代次数内无法实现最优解的探索,因此如何提高算法的逃逸能力,使其在快速定位到最优解区域同时提高最优解的精度值得我们进一步研究.鉴于此,本文提出一种多尺度协同变异的微粒群算法(MAEPSO).该算法利用不同大小方差的高斯变异机制实现解空间的探索,这种多个或大或小的变异机制,能促使整个种群以尽量分散的变异尺度来对解空间进行更加详尽的探索.同时高斯变异的范围随着适应值的变小也逐渐降低,在算法后期有利于提高最优解的精度;在利用不同大小变异算子提高局部精确搜索能力的同时,该算法同样利用均匀算子维护种群多样性.利用试验对不同评测函数进行测试均验证新算法优良的优化性能.2.解决思路提出一种多尺度协同变异的微粒群优化算法,其中算法的逃逸能力取决于不同尺度方差的高斯变异算子,不同尺度的变异有助于算法在搜索空间中进行分散式的搜索,同时变异尺度随着适应度的提升而逐渐减少,这样可以在保证逃逸能力的同时, 在算法后期提高了最优解的精度保证了算法的收敛性能.3. 实验平台本实验采用的是windows10 + Visual Studio 2013运用的win32的运行平台。

二、算法分析1.基本理论和方法我们把一个优化问题看作是在空中觅食的鸟群,那么在空中飞行的一只觅食的“鸟”就是PSO 算法在解空间中进行搜索的一个“粒子”(Particle),也是优化问题的一个解.“食物"就是优化问题的最优解.粒子在搜索空间中以一定的速度飞行,这个速度根据它本身的飞行经验和同伴的飞行经验来动态调整.所有的粒子都知道自己到目前为止发现的最好位置T，这个可以看作是粒子自己的飞行经验.此外,每个粒子还知道到目前为止整个群体中所有粒子发现的最好位置即g P (也称为gbest，gbest 是在pbest 中的最好值),这个可以看作是粒子的同伴的经验.在由这样一群随机初始化形成的粒子而组成的一个种群中,以迭代的方式进行的速度与位置的更新公式为:作为问题近似最优解。

其中, 1 r 和2 r 为两个在[0,1]范围内服从均匀分布的随机变量,用来保持群体的多样性. 1 c 和2 c 为学习因子,也称加速因子,其使粒子具有自我总结和向群体中优秀个体学习的能力,从而向自己的历史最优点以及群体的历史最优点靠近. w为惯性权重,它起着权衡局部最优能力和全局最优能力的作用.公式中右项中的三个部分分别称为：惯性部分、认知部分（代表粒子对自身的学习）、社会部分（代表着粒子间的协作）.这三个部分之间的相互平衡决定了算法的主要性能.2.具体算法通过多尺度高斯变异算子就可以保证无论最优解存在何处,算法进化到何处,都能够准确定位到最优解区域.算法具体描述如下:设尺度个数为M ,首先初始化多尺度高斯变异算子的初始方差。

初始时,方差一般设定为优化变量的取值范围,随着迭代次数的增加,多尺度高斯变异算子的方差会随之进行调整,具体调整方式如下所示,首先根据适应值的大小对种群中的微粒进行由小到大排序,然后对其进行组合,生成M 个子群，每一个子群的微粒个数为 P = N /M , K 是当前迭代次数,计算每一个子群的适应度:其中f 是微粒的适应值计算函数,不同变异尺度之间相互竞争,根据适应能力的不同而设置不同的变异能力,因此第m 个变异算子的标准差为:由于变异算子的进化是一个递归过程,排在后面的变异算子可能很大,因此对变异算子的标准差做如下规范:如果σi(k)>W / 4,则W 为待优化变量空间的宽度,重复使用上式,直到满足σi(k)<W / 4, 为了能最大范围的实现空间勘探能力,在进行完多尺度高斯变异后,同时进行一次均匀变异操作,比较后取所有变异后适应值最好的位置作为新的逃逸点,均匀变异算子实现同式(6).不同尺度算子随迭代次数的变化如图1 所示,通过本文的多尺度变异算子能实现整个搜索空间的覆盖,其中大尺度有利于实现解空间的粗搜索,可以快速定位到最优解区域,小尺度能在进化后期实现局部精确解的搜索.我们知道任何一种进化算法都存在勘探和开采两类不同的操作,PSO 算法随着迭代次数的增加粒子群会收敛于某个最优个体,粒子群个体的多样性逐渐丧失.为了防止PSO 算法陷入局部极小点,上述算法通过增加变异操作改善种群多样性来提高算法的全局解搜索能力,然而在算法的后期,变异操作往往会使微粒群不能进行精确的局部搜索无法保证算法的收敛性能.因此如何能协调勘探和开采能力是进化算法性能提高的关键.为此,本文算法采用两种方式实现勘探和开采能力间的均衡.首先,本文算法采用HPSO 的公式(4)实现PSO 进化,该公式消除了速度的惯性部分,这样就使得PSO 进化只负责整个算法的开采能力,当速度小于一定阈值时,给予速度一个变异操作——逃逸运动,描述如下:其中, f 是适应值函数,由于PSO 算法的进化公式中消除的惯性部分,使得微粒的速度快速下降,这样在有限的迭代次数内能够进行多次的逃逸运动,提高了算法的种群多样性,增强了算法全局解搜索性能.另外,由于本文算法的逃逸运动采用了多尺度高斯变异方式,较小的尺度同样能保证精确解的局部搜索能力,避免了HPSO 算法中由于下降过快且采用的均匀变异算法,使其在进化后期无法进行详细局部搜索的不足.另外,式(14)中阈值的大小将会对算法的执行效果有一定影响:如果阈值过大会导致逃逸次数增加从而打乱PSO 种群进化的原有结构,使算法无法进行局部深度搜索,开采能力有所下降;阈值过小会将导致微粒速度需要较长时间达到逃逸条件,在有限次的迭代次数内,无法进行有效地全局搜索,降低了算法搜索速度,勘探能力下降.为此结合文献[18]本文给出一个自适应的阈值设定方法.设微粒各维速度间相互独立,并对每维速度均给予一个不同的阈值,当该维度中有多个微粒速度达到这个值时,阈值就自动下降,通过自动调整每一维速度的阈值的方式达到动态调整微粒速度的目的.具体描述如下:其中, k1、k2 为常数, Size为种群大小,频率G (t) d 表示种群第d 维速度曾经发生逃逸的总次数; k1 标记调整阈值d T 和频率G (t) d 的条件值; 2 k 用来控制阈值速度下降的幅度.通过控制微粒飞行速度的方式,不仅能保证算法的全局解空间勘探性能,防止算法陷入局部最优,同时在算法后期能使粒子群进行精确的局部搜索,保证算法的收敛速度.3.算法流程图三、实验结果1.输入没有需要输入的数据，但是有三个评价函数。

如下：2.输出输出包含如下几个内容：每次更新迭代的sigma，以及最优解。

TABLET：Quadric：Rosenbrock：四、结果分析与比较1.结论本文提出一种新的多尺度协同变异的微粒群算法,利用不同大小初始方差的高斯变异机制实现解空间的探索,这种多个或大或小的变异机制,能促使整个种群以尽量分散的变异尺度来对解空间进行更加详尽的探索.同时高斯变异的范围随着适应值的变小也逐渐降低,有利于提高最优解的精度;在利用不同大小变异算子提高局部精确搜索能力的同时,该算法同样利用均匀变异算子维护种群多样性.通过6 个benchmark 数据进行测试,实验结果表明新算法能够在算法初期就能够快速定位到搜索空间的最优解区域,进而使得微粒通过进化PSO 及小尺度变异算子向着最优精确解空间逼近,使其在进化过程中保持局部最优解空间开采的能力,加快了算法的收敛速度.其中不同尺度个数对本文算法性能的影响将是本课题下一阶段研究的重点。

// MAEPSO.cpp : 定义控制台应用程序的入口点。

//#include"stdafx.h"#include<iostream>#include<math.h>#include<time.h>#include<random>#include<string.h>using namespace std;const int groupSize = 20; //种群数量const int dimNum = 30; //变量维度const int M = 5; //子群数量const int N = 6000; //迭代次数const int P = groupSize / M;const double W = 50;const int c1 = 1.4;const int c2 = 1.4;double group[groupSize][dimNum];// 粒子群double pid[groupSize][dimNum];//存放最适宜的值double pgd[dimNum];double v[groupSize][dimNum]; //粒子速度double T[dimNum]; // 每个维度上的阈值double sigma[M]; //方差int TC[dimNum]; //论文中的Gdouble w; //惯性double Tablet(double a[dimNum]) {int i = 1;double sum;sum = pow(10, 6)*a[0] * a[0];while (i < dimNum) {sum += a[i] * a[i];i++;}return sum;}double Quadric(double a[dimNum]) { //Tablet函数double sum, value = 0;for (int i = 0; i < dimNum; i++) {sum = 0;for (int j = 0; j < i; j++) {sum += a[j];}value += sum*sum;}return value;}double Rosenbrock(double a[dimNum]) {double value = 0;for (int i = 0; i < dimNum - 1; i++) {value += 100 * (a[i + 1] - a[i] * a[i])*(a[i + 1] - a[i] * a[i]) + (a[i] - 1)*(a[i] - 1);}return value;}double FitX(int m) { //求子群的平均适应度double sum = 0;for (int i = 0; i < P; i++) {sum += Tablet(group[m*P + i]);}sum = sum / P;return sum;}double MaxFitX() { //求所有子群的平均适应度的最大值double max = DBL_MIN;for (int i = 0; i < M; i++) {if (max < FitX(i)) {max = FitX(i);}}return max;}double MinFitX() { //求所有子群的平均适应度的最小值double min = DBL_MAX;for (int i = 0; i < M; i++) {if (min > FitX(i)) {min = FitX(i);}}return min;}void Gsort() { //对粒子进行排序double temp[dimNum];for (int i = 0; i < groupSize; i++) {for (int j = i; j < groupSize; j++) {if (Tablet(group[i]) > Tablet(group[j])) {memcpy(temp, group[i], sizeof(temp));memcpy(group[i], group[j], sizeof(group[i]));memcpy(group[j], temp, sizeof(temp));memcpy(temp, pid[i], sizeof(temp));memcpy(pid[i], pid[j], sizeof(pid[i]));memcpy(pid[j], temp, sizeof(temp));}}}}void updateSigma() { //论文中的2.1double sum = 0;for (int i = 0; i < M; i++) {sum += FitX(i);}for (int i = 0; i < M; i++) {sigma[i] = sigma[i] * exp((M*FitX(i) - sum) / (MaxFitX() - MinFitX()));while (sigma[i] > W / 4) {sigma[i] = abs(W / 4 - sigma[i]);}cout << sigma[i] <<" ";}}double sampleNormal() { //随机符合N[0,1]的正态分布double u = ((double)rand() / (RAND_MAX)) * 2 - 1;double v = ((double)rand() / (RAND_MAX)) * 2 - 1;double r = u * u + v * v;if (r == 0 || r > 1) return sampleNormal();double c = sqrt(-2 * log(r) / r);return u * c;}void updateSpeed() { //更新速度V和粒子群double Vmax, minf, mintemp, temp;for (int i = 0; i < M; i++) { //迭代每一个粒子群for (int j = 0; j < P; j++) { //迭代每一个粒子for (int z = 0; z < dimNum; z++) { //迭代每一个维度if (v[i*P + j][z] == 0) { //HPSO 中的(4)double r1 = rand()*1.0 / RAND_MAX;int r2 = 1;if (r1 > 0.5) {r2 = -1;}v[i*P + j][z] = r1*r2*pow(0.1, 12);}if (v[i*P + j][z] < T[z]) { //论文中的（14）TC[z]++;temp = group[i*P + j][z]; //用来暂存更新前的Xidminf = DBL_MAX; //最小适应度for (int k = 0; k < M; k++) {group[i*P + j][z] = temp + sampleNormal()*sigma[k];if (minf > Tablet(group[i*P + j])) {minf = Tablet(group[i*P + j]);mintemp = group[i*P + j][z];}}Vmax = W - abs(temp);group[i*P + j][z] = temp + (rand()*2.0 / RAND_MAX - 1)*Vmax;if (minf < Tablet(group[i*P + j])) {v[i*P + j][z] = mintemp - temp;}else {v[i*P + j][z] = group[i*P + j][z] - temp;}group[i*P + j][z] = temp;}v[i*P + j][z] = w*v[i*P + j][z] + c1*(rand()*1.0 / RAND_MAX)*(pid[i*P + j][z] - group[i*P + j][z]) + c2*(rand()*1.0 / RAND_MAX)*(pgd[z] - group[i*P + j][z]); //PSO group[i*P + j][z] += v[i*P + j][z]; //进行种群更新if (Tablet(group[i*P + j]) < Tablet(pid[i*P + j])) {memcpy(pid[i*P + j], group[i*P + j], sizeof(group[i*P + j]));if (Tablet(pid[i*P + j]) < Tablet(pgd))memcpy(pgd, pid[i*P + j], sizeof(pgd));}}}}}void updateT() {int k1 = 5, k2 = 10;for (int i = 0; i < dimNum; i++) {if (TC[i] > k1){T[i] = T[i] / k2;TC[i] = 0;}}}void initial() {for (int i = 0; i < groupSize; i++){for (int j = 0; j < dimNum; j++) {group[i][j] = (rand()*2.0 / RAND_MAX - 1)*W;pid[i][j] = group[i][j];v[i][j] = 0;}}Gsort();memcpy(pgd, pid[0], sizeof(pgd));for (int i = 0; i < M; i++) {sigma[i] = 2*W;}for (int i = 0; i < dimNum; i++) {T[dimNum] = 0.5;TC[dimNum] = 0;}}int main(){srand(time(0));initial();for (int i = 0; i < N; i++) {w = 0.95 - 0.55*i / N;cout <<"第"<< i <<"次迭代的SIGMA为：";updateSigma();updateSpeed();updateT();cout <<"得出最优解适应度为：";cout << Tablet(pgd) << endl;}cout <<"得出最优解为："<< endl; for (int j = 0; j < dimNum; j++) { cout << pgd[j] << endl;}cout <<"得出最优解适应度为：";cout << Tablet(pgd) << endl;return 0;}。