用MATLAB编写PSO算法及实例

1、用 MATLAB编写 PSO算法及实例用 MATLAB编写 PSO算法及实例1.1 粒子群算法PSO从这种模型中得到启示并用于解决优化问题。 PSO 中，每个优化问题的 潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化 的函数决定的适值 ( fitness value) ，每个粒子还有一个速度决定它们飞翔的 方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。PSO初始化为一群随机粒子 (随机解 ) ，然后通过迭代找到最优解。 在每一次 迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最 优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解

2、，这个 极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻 居，那么在所有邻居中的极值就是局部极值。假设在一个 D维的目标搜索空间中， 有N个粒子组成一个群落， 其中第i个 粒子表示为一个 D 维的向量X i (xi1,xi2, ,xiD ) ， i 1,2, ,N 。 ，。第i 个粒子的“飞行 ”速度也是一个 D 维的向量，记为Vi (vi1,vi2， ,viD ) ， i 1,2, 3 。 ，。第i 个粒子迄今为止搜索到的最优位置称为个体极值，记为pbest (pi1,pi2, ,piD)，i 1,2, ,N 。整个粒子群迄今为止搜索到的最优位置为全局极值，记为gbest

3、 (pg1, pg2, , pgD )在找到这两个最优值时，粒子根据如下的公式 (1.1) 和( 1.2) 来更新自己的速度和位置：(1.1)vid w vid c1r1 pid xidc2 r2 ( pgd xid )(1. 2)xidxidvid其中： c1和c2为学习因子，也称加速常数 (acceleration constant)，r1和r2为0，1 范围内的均匀随机数。式 (1.1) 右边由三部分组成，第一部分为“惯性 (inertia) ”或“动量 (momentum)”部分，反映了粒子的运动“习惯 (habit) ”， 代表粒子有维持自己先前速度的趋势；第二部分为“认知 (cog

4、nition) ”部分， 反映了粒子对自身历史经验的记忆 (memory)或回忆 (remembrance) ，代表粒子有 向自身历史最佳位置逼近的趋势；第三部分为“社会 (social) ”部分，反映了 粒子间协同合作与知识共享的群体历史经验。二、算法设计2.1 算法流程图2.2 算法实现算法的流程如下： 初始化粒子群，包括群体规模 N ，每个粒子的位置 xi和速度 Vi 计算每个粒子的适应度值 Fit i；%它的适应度就是指 目标 函数的值 。一般来说，目标函数的选择由具体问 题来决定， 假如是背包问题，适应度即放入包中 物体的总价格。 初始粒子位置和速度的位置一 般随机产生。 但是在某些

5、领域，如果已有其他的 算法可以产生可行解的话， 可以用这个可行解来 初始化，这样更容易得到最优的解 对每个粒子，用它的适应度值 Fiti 和个体极值 pbe（st i）比较，如果Fiti pbest（i） ，则用 Fiti 替换掉 pbe（st i）； 对每个粒子，用它的适应度值 Fiti和全局极值gbest比较，如果Fiti pbest （i ）则用 Fiti替gbest； 根据公式（ 1.1 ），（ 1.2 ）更新粒子的速度 vi 和位置 xi ； 如果满足结束条件 （ 误差足够好或到达最大循环次数 ） 退出，否则返回。2.3 参数选择本算法中主要的参数变量为 w（惯性权值），c1 ，c2

6、（加速因子），N （种群数），M （迭代次数）， D （粒子维数）1）种群规模通常，种群太小则不能提供足够的采样点，以致算法性能很差；种群太大 尽管可以增加优化信息，阻止早熟收敛的发生，但无疑会增加计算量，造成收 敛时间太长，表现为收敛速度缓慢。种群规模一般设为1001000。本文选择种群规模为 100。（2）最大迭代次数 迭代次数越多能保证解的收敛性，但是影响运算速度，本文选1000次。（3）惯性权值惯性权重 w 表示在多大程度上保留原来的速度。 w较大，全局收敛能力强， 局部收敛能力弱； w较小，局部收敛能力强，全局收敛能力弱。本文选 0.6 。（4）加速因子加速常数 c2和c2 分别用于

7、控制粒子指向自身或邻域最佳位置的运动。文献20 建议 c1 c2 4.0 ，并通常取 c1 c2 2 。本文也取 c1 c2 2 。（5）粒子维数本文中粒子维数取决于待优化函数的维数。需要说明的是，本文的程序允许改变这些参数，因为本文编写的程序参照 matlab 工具箱，留给用户解决这类问题一个接口函数，上述的各个参数正是接 口函数的参数，因此允许改变。另外对于 w 和c也可采用变参数法，即随迭代次 数增加，利用经验公式使它们动态调整，本文采用固定值。3.1 求三维函数 f=x（1）.2+x（2）.2+x（3）.2 的最小值步骤： 1. 初始化 x,v; 2. 求出每个粒子的适应 值；3. 初

8、始化 pb,pg 个体最优和全局最优； 4. 根据式子更新 x,v; 5. 是否满足条件，满足跳出循环，否则重复 2-4 步尝试编码：1)pso.m 文件%此算法是 PSO 算法，汪汪的 20161024 号版本function xm,fv=PSO(fitness,N,c1,c2,w,M,D)%xm,fv 算法最后得到的最优解时的 x 及最优解， fitness 为适应度，即要优化的目标函数，N为种群数量， c1,c2 为学习因子， w 为惯性权重， M为迭代次数， D 为粒子的维数%format long ;%初始化种群 for i=1:Nfor j=1:Dx(i,j)=randn;%随机初

9、始化位置v(i,j)=randn;%随机初始化速度endend%先计算各个粒子的适应度pi ，并初始化 y- 粒子个体极值， pg- 全局极值for i=1:Np(i)=fitness(x(i,:);%适应度 问题：将 x(i,:) 改成 x(i,j) 是否可以，答不能y(i,:)=x(i,:);%个体极值endpg=x(N,:);%初始化全局极值 / 最优for i=1:N-1if fitness(x(i,:)fitness(pg)pg=x(i,:);%替换并选出全局极值endend%进入粒子群算法主要循环，更新v 及 xfor t=1:Mfor i=1:Nv(i,:)=w*v(i,:)+c

10、1*rand*(y(1,:)-x(1,:)+c2*rand*(pg-x(i,:);x(i,:)=x(i,:)+v(i,:);if fitness(x(i,:)p(i)p(i)=fitness(x(i,:);y(i,:)=x(i,:);endif p(i)fitness(pg)pg=y(i,:);end%M次迭代后最优解endpbest(t)=fitness(pg);endxm=pg; %为何要共轭转置？ fv=fitness(pg);(2)目标函数 fitness.m 文件function f=fitness(x) f=x(1).2+x(2).2+x(3).2 ;end需要说明的是，针对不同的函数优化，只需要改变目标函数就可以。 (3)在命令行输入或建立调用 m文件在命令行先后输入 xm,fv = PSO(fitness,100,2,2,0.6,1000,3) 立包涵该语句的 m文件，运行即可得到结果。四、 结果与分析xm=1.0e-04 *-0.2857301402295650.676783696397148-0.250529540096653 fv =6.024429352056337e-09fv 是最优值， xm为最优值对应的自变量值。3.2 高斯函数，或建%生x y=me