《粒子群算法vba》doc版.docx

水资源系统分析理论与应用课程论文学 院： 水文水资源 专 业： 水利工程 题 目： 粒子群算法（VBA） 学 号： 171601010078 指导老师： 王建群 姓 名： 张轶文 2017年11月14日粒子群算法张轶文（学号：161601010078，专业：水利工程）一、粒子群算法的起源粒子群优化(Particle Swarm Optimization, PSO)算法是Kennedy和Eberhart受人工生命研究结果的启发、通过模拟鸟群觅食过程中的迁徙和群聚行为而提出的一种基于群体智能的全局随机搜索算法，1995年IEEE国际神经网络学术会议发表了题为“Particle Swarm Optimization”的论文，标志PSO算法诞生(注：国内也有很多学者译为“微粒群优化”)。它与其他进化算法一样，也是基于“种群”和“进化”的概念，通过个体间的协作与竞争，实现复杂空间最优解的搜索；同时，PSO又不像其他进化算法那样对个体进行交叉、变异、选择等进化算子操作，而是将群体(swarm)中的个体看作是在D维搜索空间中没有质量和体积的粒子(particle)，每个粒子以一定的速度在解空间运动，并向自身历史最佳位置pbest和邻域历史最佳位置pbest聚集，实现对候选解的进化。PSO算法具有很好的生物社会背景而易理解、参数少而易实现，对非线性、多峰问题均具有较强的全局搜索能力，在科学研究与工程实践中得到了广泛关注。 自然界中各种生物体均具有一定的群体行为，而人工生命的主要研究领域之一是探索自然界生物的群体行为，从而在计算机上构建其群体模型。自然界中的鸟群和鱼群的群体行为一直是科学家的研究兴趣，生物学家Craig Reynolds在1987年提出了一个非常有影响的鸟群聚集模型，在他的仿真中，每一个个体遵循： (1) 避免与邻域个体相冲撞； (2) 匹配邻域个体的速度； (3) 飞向鸟群中心，且整个群体飞向目标。 仿真中仅利用上面三条简单的规则，就可以非常接近的模拟出鸟群飞行的现象。1990年，生物学家Frank Heppner也提出了鸟类模型，它的不同之处在于：鸟类被吸引飞到栖息地。在仿真中，一开始每一只鸟都没有特定的飞行目标，只是使用简单的规则确定自己的飞行方向和飞行速度（每一只鸟都试图留在鸟群中而又不相互碰撞），当有一只鸟飞到栖息地时，它周围的鸟也会跟着飞向栖息地，这样，整个鸟群都会落在栖息地。 1995年，美国社会心理学家James Kennedy和电气工程Russell Eberhart共同提出了粒子群算法，其基本思想是受对鸟类群体行为进行建模与仿真的研究结果的启发。他们的模型和仿真算法主要对Frank Heppner的模型进行了修正，以使粒子飞向解空间并在最好解处降落。Kennedy在他的书中描述了粒子群算法思想的起源。自20世纪30年代以来，社会心理学的发展揭示：我们都是鱼群或鸟群聚集行为的遵循者。在人们的不断交互过程中，由于相互的影响和模仿，他们总会变得更相似，结果就形成了规范和文明。人类的自然行为和鱼群及鸟群并不类似，而人类在高维认知空间中的思维轨迹却与之非常类似。思维背后的社会现象远比鱼群和鸟群聚集过程中的优美动作复杂的多：首先，思维发生在信念空间，其维数远远高于3；其次，当两种思想在认知空间会聚于同一点时，我们称其一致，而不是发生冲突。二、算法的基本原理 PSO 中，每个优化问题的潜在解都是搜索空间中的一只鸟，称之为粒子。所有的粒子都有一个由被优化的函数决定的适值(fitness value) ，每个粒子还有一个速度决定它们飞翔的方向和距离。然后粒子们就追随当前的最优粒子在解空间中搜索。PSO初始化为一群随机粒子(随机解)，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己；第一个就是粒子本身所找到的最优解，这个解称为个体极值；另一个极值是整个种群目前找到的最优解，这个极值是全局极值。另外也可以不用整个种群而只是用其中一部分作为粒子的邻居，那么在所有邻居中的极值就是局部极值。 假设在一个维的目标搜索空间中，有m个粒子组成一个群落，其中第个粒子表示为一个维的向量，第个粒子的“飞行 ”速度也是一个维的向量，记为 ，第个粒子迄今为止搜索到的最优位置称为个体极值，记为 ，整个粒子群迄今为止搜索到的最优位置为全局极值，记为 在找到这两个最优值时，粒子根据如下的公式(2.1)和( 2.2)来更新自己的速度和位置： (2.1) (2. 2)其中：和为学习因子，也称加速常数，和为0，1范围内的均匀随机数。式(2.1)右边由三部分组成，第一部分为“惯性”或“动量”部分，反映了粒子的运动“习惯”，代表粒子有维持自己先前速度的趋势；第二部分为“认知”部分，反映了粒子对自身历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；第三部分为“社会”部分，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻域历史最佳位置逼近的趋势，根据经验，通常。是粒子的速度。，是常数，由用户设定用来限制粒子的速度。和是介于0,1之间的随机数。带有惯性权重的改进粒子群算法其进化过程为： （2.3） （2.4）在式(2.1)中，第一部分表示粒子先前的速度，用于保证算法的全局收敛性能；第二部分、第三部分则是使算法具有局部收敛能力。可以看出，式（2.3)中惯性权重w表示在多大程度上保留原来的速度。w较大，全局收敛能力强，局部收敛能力弱；较小，局部收敛能力强，全局收敛能力弱。对w的取值，目前多采用线性动态惯性权重，更新公式如下：W(t)=Wmax-(Wmax-Wmin)/Tmax*t三、算法的基本步骤基本粒子群的算法步骤如下：1、Initial：初始化粒子群体（群体规模为n），包括随机位置和速度。2、Evaluation：根据fitness function，评价每个例子的适应度。3、 Find the Pbest：对每个粒子，将其当前适应值与其个体历史最佳位置（pbest）对应的适应值作比较，如果当前的适应值更高，则将用当前位置更新历史最佳位置pbest。4、Find the Gbest：对每个粒子，将其当前适应值与全局最佳位置（gbest）对应的适应值作比较，如果当前的适应值更高，则将用当前粒子的位置更新全局最佳位置gbest。5、Update the Veloctiy：根据公式更新每个粒子的速度与位置。 图1 粒子群算法的流程图四、函数介绍（1）f1：Sphere Model function这里，为函数的维数。该函数是凸函数，有唯一的全局最小值（2）f2：Goldstein-Price function函数该函数有若干个极小值点，全局最小值为:（3）f3：Schwefel function这里，为函数的维数。该函数的全局最小值为:，该函数的全局最优点在搜索空间的边缘，并与次最优解相距甚远，搜索算法极易收敛到错误的方向。（4）f4：Rana function该函数的全局最小值为:，该函数是一个具有很多极小值点的函数，最优解位于搜索空间的角落。（5）f5：Ackley function这里，为函数的维数。该函数的全局最小值为（6）f6：Rosenbrock函数该函数的全局最小值为，该函数的全局最优点隐藏在一个狭长的、抛物线形的平直峡谷里，比较难于找到。（7）f7：Rastrigin函数该函数的全局最小值为，该函数的全局最优点被许多局部极小值点所包围，非常难于找到。（8）f8：Griewank函数这里及以下，表示函数的维数。该函数的全局最小值为，该函数具有许多分布很广泛的局部极小值点，而且这些局部极小值相差不大，所以极不易收敛到全局最优解；（9）f9：多个全局最优点函数 有6个局部极小点, 全局极小点为（- 0.0898，0.7126）和(0.0898，-0.7126)，最小值为 1.031628。（10）f10：Branin funciton有3个全局最小值点：五、测试结果1、计算精度展示当取粒子个数i=100,函数维度=2，学习因子=2。结果如下:表1.1全局粒子群最优解函数理论最优值第一次计算第二次计算第三次计算f1最优值00.0020.0000.000迭代次数101001000f2最优值33.0003.0003.000迭代次数101001000f3最优值-837.9658-937.1790-828.6437-837.9658迭代次数101001000f4最优值-511.733-511.708-511.708-511.708迭代次数101001000f5最优值04.8314.6714.671迭代次数101001000f6最优值00.0150.0000.000迭代次数101001000f7最优值00.0000.0000.000迭代次数101001000f8最优值0-0.295-0.299-0.287迭代次数101001000f9最优值-1.031628-1.031599-1.031599-1.031599迭代次数101001000f10最优值0.3978870.3980990.3980000.398099迭代次数1010810001.根据粒子群算法获得上表中的结果，结果相对精确，误差较小。2.维度限制大，只做了二维，结果相对单一，无法推测多维度的结果。六、总结首先我利用了二个月的时间学习了VBA编程，Visual Basic for Applications（VBA）是Visual Basic的一种宏语言，是微软开发出来在其桌面应用程序中执行通用的自动化(OLE)任务的编程语言。主要能用来扩展Windows的应用程序功能，特别是Microsoft Office软件。初步学习VBA中的定义语句、选择语句、判断语句和循环语句，能自己编写一定难度的程序。在本次粒子群编程中，因为时间有限，所以仅仅编写了二维的粒子群算法。通过十个测试函数，不但学习认识了粒子群算法的起源、原理和基本步骤，更强化了个人VBA编程能力，提升编程熟练度。就如王建群老师说的，现在都21世纪了，情感机器人都出现了，我们怎么能不会编程。在以后的学习中，本人会更加积极主动的学习编程，学习水文知识。附件F1Option ExplicitSub Text()Dim i As Integer, j As Integer, l As Integer, w As SingleDim c1 As Integer, c2 As Integer, x1(100) As SingleDim x2(100) As Single, x3(100) As Single, v1(100) As SingleDim v2(100) As Single, v3(100) As Single, vmax As SingleDim vmin As Single, pbest1(100) As Single, pbest2(100) As SingleDim pbest3(100) As Single, gbest1 As Single, gbest2 As SingleDim gbest3 As Single, y(100) As Single, a As Single, b As SingleRange(A1).Value = X1Range(B1).Value = V1Range(C1).Value = X2Range(D1).Value = V2Range(E1).Value = X3Range(F1).Value = V3Range(G1).Value = pbest1Range(H1).Value = pbest2Range(I1).Value = pbest3Range(J1).Value = gbest1Range(K1).Value = gbest2Range(L1).Value = gbest3w = 0.4c1 = 2c2 = 2vmax = 10.24vmin = -10.24gbest1 = 100For i = 0 To 99x1(i) = Rnd() * 10.24 - 5.12v1(i) = Rnd() * 10.24Range(A & (i + 2).Value = Round(x1(i), 3)Range(B & (i + 2).Value = Round(v1(i), 3)Next igbest2 = 100For i = 0 To 99x2(i) = Rnd() * 10.24 - 5.12v2(i) = Rnd() * 10.24Range(C & (i + 2).Value = Round(x2(i), 3)Range(D & (i + 2).Value = Round(v2(i), 3)Next igbest3 = 100For i = 0 To 99x3(i) = Rnd() * 10.24 - 5.12v3(i) = Rnd() * 10.24Range(E & (i + 2).Value = Round(x3(i), 3)Range(F & (i + 2).Value = Round(v3(i), 3)Next ib = 100For i = 0 To 99pbest1(i) = Round(x1(i), 3)Range(G & (i + 2).Value = pbest1(i)a = x1(i) * x1(i)If a b Thenb = agbest1 = pbest1(i)End IfNext iRange(J2).Value = gbest1b = 100For i = 0 To 99pbest2(i) = Round(x2(i), 3)Range(H & (i + 2).Value = pbest2(i)a = x2(i) * x2(i)If a b Thenb = agbest2 = x2(i)End IfNext iRange(K2).Value = gbest2b = 100For i = 0 To 99pbest3(i) = Round(x3(i), 3)Range(I & (i + 2).Value = pbest3(i)a = x3(i) * x3(i)If a b Thenb = agbest3 = x3(i)End IfNext iRange(L2).Value = gbest3For j = 0 To 99For i = 0 To 99v1(i) = w * v1(i) + c1 * Rnd() * (pbest1(i) - x1(i) + c2 * Rnd() * (gbest1 - x1(i)x1(i) = x1(i) + v1(i)If x1(i) vmax) Thenv1(i) = vmaxEnd IfIf (v1(i) 5.12) Thenx1(i) = 5.12End IfIf (x1(i) -5.12) Thenx1(i) = -5.12End IfNext iFor i = 0 To 99v2(i) = w * v2(i) + c1 * Rnd() * (pbest2(i) - x2(i) + c2 * Rnd() * (gbest2 - x2(i)x2(i) = x2(i) + v2(i)If x2(i) vmax) Thenv2(i) = vmaxEnd IfIf (v2(i) 5.12) Thenx2(i) = 5.12End IfIf (x2(i) -5.12) Thenx2(i) = -5.12End IfNext iFor i = 0 To 99v3(i) = w * v3(i) + c1 * Rnd() * (pbest3(i) - x3(i) + c2 * Rnd() * (gbest3 - x3(i)x3(i) = x3(i) + v3(i)If x3(i) vmax) Thenv3(i) = vmaxEnd IfIf (v3(i) 5.12) Thenx3(i) = 5.12End IfIf (x3(i) -5.12) Thenx3(i) = -5.12End IfNext iFor i = 0 To 99If x1(i) * x1(i) pbest1(i) * pbest1(i) Thenpbest1(i) = x1(i)End IfIf pbest1(i) * pbest1(i) gbest1 * gbest1 Thengbest1 = pbest1(i)End IfIf x2(i) * x2(i) pbest2(i) * pbest2(i) Thenpbest2(i) = x2(i)End IfIf pbest2(i) * pbest2(i) gbest2 * gbest2 Thengbest2 = pbest2(i)End IfIf x3(i) * x3(i) pbest3(i) * pbest3(i) Thenpbest3(i) = x3(i)End IfIf pbest3(i) * pbest3(i) gbest3 * gbest3 Thengbest3 = pbest3(i)End IfNext iNext ja = 100For i = 0 To 99y(i) = x1(i) * x1(i) + x2(i) * x2(i) + x3(i) * x3(i)If y(i) a Thena = Round(y(i), 3)End IfNext iRange(M1).Value = X1Range(N1).Value = X2Range(O1).Value = X3Range(P1).Value = YRange(Q1).Value = YminRange(Q2).Value = aFor i = 0 To 99Range(M & (i + 2).Value = Round(x1(i), 3)Range(N & (i + 2).Value = Round(x2(i), 3)Range(O & (i + 2).Value = Round(x3(i), 3)Range(P & (i + 2).Value = Round(y(i), 3)Next iEnd SubF2Option ExplicitSub Text()Dim i As IntegerDim j As IntegerDim w As SingleDim c1 As IntegerDim c2 As IntegerDim x1(150) As SingleDim x2(150) As SingleDim v1(150) As SingleDim v2(150) As SingleDim y(150) As SingleDim vmax As SingleDim vmin As SingleDim pbest1(150) As SingleDim pbest2(150) As SingleDim gbest1 As SingleDim gbest2 As SingleDim a As SingleDim y1 As SingleDim y2 As SingleDim y3 As SingleDim y4 As SingleRange(A1).Value = X1Range(B1).Value = V1Range(C1).Value = X2Range(D1).Value = V2Range(E1).Value = pbest1Range(F1).Value = pbest2Range(G1).Value = gbest1Range(H1).Value = gbest2Range(I1).Value = Yw = 0.01c1 = 2c2 = 2a = 1000vmax = 4vmin = -4gbest1 = 2gbest2 = 2For i = 0 To 99x1(i) = Rnd() * 4 - 2v1(i) = Rnd() * 4x2(i) = Rnd() * 4 - 2v2(i) = Rnd() * 4Range(A & (i + 2).Value = Round(x1(i), 3)Range(B & (i + 2).Value = Round(v1(i), 3)Range(C & (i + 2).Value = Round(x2(i), 3)Range(D & (i + 2).Value = Round(v2(i), 3)Next iFor i = 0 To 99y1 = x1(i) + x2(i) + 1y2 = 19 - 14 * x1(i) + 3 * x1(i) * x1(i) - 14 * x2(i) + 6 * x1(i) * x2(i) + 3 * x2(i) * x2(i)y3 = 2 * x1(i) - 3 * x2(i)y4 = 18 - 32 * x1(i) + 12 * x1(i) * x1(i) + 48 * x2(i) - 36 * x1(i) * x2(i) + 27 * x2(i) * x2(i)y(i) = (1 + y1 * y1 * y2) * (30 + y3 * y3 * y4)Range(I & (i + 2).Value = Round(y(i), 3)pbest1(i) = Round(x1(i), 3)pbest2(i) = Round(x2(i), 3)Range(E & (i + 2).Value = pbest1(i)Range(F & (i + 2).Value = pbest2(i)If y(i) a Thena = y(i)gbest1 = x1(i)gbest2 = x2(i)End IfNext iRange(G2).Value = gbest1Range(H2).Value = gbest2For j = 0 To 9For i = 0 To 99v1(i) = w * v1(i) + c1 * Rnd() * (pbest1(i) - x1(i) + c2 * Rnd() * (gbest1 - x1(i)x1(i) = x1(i) + v1(i)v2(i) = w * v2(i) + c1 * Rnd() * (pbest2(i) - x2(i) + c2 * Rnd() * (gbest2 - x2(i)x2(i) = x2(i) + v2(i)y1 = x1(i) + x2(i) + 1y2 = 19 - 14 * x1(i) + 3 * x1(i) * x1(i) - 14 * x2(i) + 6 * x1(i) * x2(i) + 3 * x2(i) * x2(i)y3 = 2 * x1(i) - 3 * x2(i)y4 = 18 - 32 * x1(i) + 12 * x1(i) * x1(i) + 48 * x2(i) - 36 * x1(i) * x2(i) + 27 * x2(i) * x2(i)If (1 + y1 * y1 * y2) * (30 + y3 * y3 * y4) y(i) Thenpbest1(i) = x1(i)pbest2(i) = x2(i)y(i) = (1 + y1 * y1 * y2) * (30 + y3 * y3 * y4)End IfIf y(i) a Thena = y(i)gbest1 = pbest1(i)gbest2 = pbest2(i)End IfNext iNext ja = 1000For i = 0 To 99y1 = x1(i) + x2(i) + 1y2 = 19 - 14 * x1(i) + 3 * x1(i) * x1(i) - 14 * x2(i) + 6 * x1(i) * x2(i) + 3 * x2(i) * x2(i)y3 = 2 * x1(i) - 3 * x2(i)y4 = 18 - 32 * x1(i) + 12 * x1(i) * x1(i) + 48 * x2(i) - 36 * x1(i) * x2(i) + 27 * x2(i) * x2(i)y(i) = (1 + y1 * y1 * y2) * (30 + y3 * y3 * y4)Range(J & i + 2).Value = Round(x1(i), 3)Range(K & i + 2).Value = Round(x2(i), 3)Range(L & i + 2).Value = Round(y(i), 3)If y(i) a Thena = Round(y(i), 5)End IfNext iRange(M1).Value = gbestRange(N1).Value = gbestRange(M2).Value = gbest1Range(N2).Value = gbest2Range(J1).Value = X1Range(K1).Value = X2Range(L1).Value = YRange(O1).Value = YminRange(O2).Value = aEnd SubF3Option ExplicitSub Text()Dim i As IntegerDim j As IntegerDim w As SingleDim c1 As IntegerDim c2 As IntegerDim x1(100) As SingleDim x2(100) As SingleDim v1(100) As SingleDim v2(100) As SingleDim y(100) As SingleDim vmax As SingleDim vmin As SingleDim pbest1(100) As SingleDim pbest2(100) As SingleDim gbest1 As SingleDim gbest2 As SingleDim a As SingleRange(A1).Value = X1Range(B1).Value = V1Range(C1).Value = X2Range(D1).Value = V2Range(E1).Value = pbest1Range(F1).Value = pbest2Range(G1).Value = gbest1Range(H1).Value = gbest2Range(I1).Value = Yw = 0.4c1 = 2c2 = 2a = 1000vmax = 1000vmin = -1000gbest1 = 500gbest2 = 500For i = 0 To 99x1(i) = Rnd() * 1000 - 500v1(i) = Rnd() * 1000x2(i) = Rnd() * 1000 - 500v2(i) = Rnd() * 1000Range(A & (i + 2).Value = Round(x1(i), 3)Range(B & (i + 2).Value = Round(v1(i), 3)Range(C & (i + 2).Value = Round(x2(i), 3)Range(D & (i + 2).Value = Round(v2(i), 3)Next iFor i = 0 To 99y(i) = -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i)Range(I & (i + 2).Value = Round(y(i), 3)pbest1(i) = Round(x1(i), 3)pbest2(i) = Round(x2(i), 3)Range(E & (i + 2).Value = pbest1(i)Range(F & (i + 2).Value = pbest2(i)If y(i) vmax) Thenv1(i) = vmaxEnd IfIf (v1(i) 500) Thenx1(i) = 500End IfIf (x1(i) -500) Thenx1(i) = -500End IfIf -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i) y(i) Thenpbest1(i) = x1(i)y(i) = -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i)End IfIf y(i) vmax) Thenv2(i) = vmaxEnd IfIf (v2(i) 500) Thenx2(i) = 500End IfIf (x2(i) -500) Thenx2(i) = -500End IfIf -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i) Thenpbest2(i) = x2(i)y(i) = -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i)End IfIf y(i) a Thena = y(i)gbest2 = pbest2(i)End IfNext iNext ja = 1000For i = 0 To 99y(i) = -x1(i) * Sin(Sqr(Abs(x1(i) - x2(i) * Sin(Sqr(Abs(x2(i)Range(J & i + 2).Value = Round(x1(i), 3)Range(K & i + 2).Value = Round(x2(i), 3)Range(L & i + 2).Value = Round(y(i), 3)If y(i) a Thena = Round(y(i), 4)End IfNext iRange(M1).Value = gbestRange(N1).Value = gbestRange(M2).Value = gbest1Range(N2).Value = gbest2Range(J1).Value = X1Range(K1).Value = X2Range(L1).Value = YRange(O1).Value = YminRange(O2).Value = aEnd SubF4Option ExplicitSub Text()Dim i As IntegerDim j As IntegerDim w As SingleDim c1 As IntegerDim c2 As IntegerDim x1(150) As SingleDim x2(150) As SingleDim v1(150) As SingleDim v2(150) As SingleDim y(150) As SingleDim vmax As SingleDim vmin As SingleDim pbest1(150) As SingleDim pbest2(150) As SingleDim gbest1 As SingleDim gbest2 As SingleDim a As SingleDim y1 As SingleDim y2 As SingleRange(A1).Value = X1Range(B1).Value = V1Range(C1).Value = X2Range(D1).Value = V2Range(E1).Value = pbest1Range(F1).Value = pbest2Range(G1).Value = gbest1Range(H1).Value = gbest2Range(I1).Value = Yw = 1.4c1 = 2c2 = 2a = 1000vmax = 1024vmin = -1024gbest1 = 512gbest2 = 512For i = 0 To 149x1(i) = Rnd() * 1024 - 512v1(i) = Rnd() * 1024x2(i) = Rnd() * 1024 - 512v2(i) = Rnd() * 1024Range(A & (i + 2).Value = Round(x1(i), 3)Range(B & (i + 2).Value = Round(v1(i), 3)Range(C & (i + 2).Value = Round(x2(i), 3)Range(D & (i + 2).Value = Round(v2(i), 3)Next iFor i = 0 To 99y1 = x1(i) * Sin(Sqr(Abs(x2(i) + 1 - x1(i) * Cos(Sqr(Abs(x2(i) + 1 + x1(i)y2 = (x2(i) + 1) * Cos(Sqr(Abs(x2(i) + 1 - x1(i) * Sin(Sqr(Abs(x2(i) + 1 + x1(i)y(i) = y1 + y2Range(I & (i + 2).Value = Round(y(i), 3)pbest1(i) = Round(x1(i), 3)pbest2(i) = Round(x2(i), 3)Range(E & (i + 2).Value = pbest1(i)Range(F & (i + 2).Value = pbest2(i)If y(i) a Thena = y(i)gbest1 = x1(i)gbest2 = x2(i)End IfNext iRange(G2).Value = gbest1Range(H2).Value = gbest2For j = 0 To 9For i = 0 To 99v1(i) = w * v1(i) + c1 * Rnd() * (pbest1(i) - x1(i) + c2 * Rnd() * (gbest1 - x1(i)x1(i) = x1(i) + v1(i)v2(i) = w * v2(i) + c1 * Rnd() * (pbest2(i) - x2(i) + c2 * Rnd() * (gbest2 - x2(i)x2(i) = x2(i) + v2(i)y1 = x1(i) * Sin(Sqr(Abs(x2(i) + 1 - x1(i) * Cos(Sqr(Abs(x2(i) + 1 + x1(i)y2 = (x2(i) + 1) * Cos(Sqr(Abs(x2(i) + 1 - x1(i) * Sin(Sqr(Abs(x2(i) + 1 + x1(i)If y1 + y2 vmax) Thenv1(i) = vmaxEnd IfIf (v1(i) 512) Thenx1(i) = 512End IfIf (x1(i) vmax) Thenv2(i) = vmaxEnd IfIf (v2(i) 512) Thenx2(i) = 512End IfIf (x2(i) -512) Thenx2(i) = -512End IfIf y(i) a Thena = y(i)gbest1 = pbest1(i)gbest2 = pbest2(i)End IfNext iNext ja = 1000For i = 0 To 99y1 = x1(i) * Sin(Sqr(Abs(x2(i) + 1 - x1(i) * Cos(Sqr(Abs(x2(i) + 1 + x1(i)y2 = (x2(i) + 1) * Cos(Sqr(Abs(x2(i) + 1 - x1(i) * Sin(Sqr(Abs(x2(i) + 1 + x1(i)y(