中规中矩, 没有ABC快
1 -- 粒子群算法 PSO(全局,变权重) 2 3 -- lua 相关简易操作 4 sin = math.sin 5 cos = math.cos 6 sqrt = math.sqrt 7 pi = math.pi 8 random = math.random 9 exp = math.exp 10 int = math.floor 11 12 13 -- 获得不同的随机序列 14 math.randomseed(os.time()) 15 16 17 -- ==============================求解函数============================== 18 dim = 2 -- 解的维数 19 20 domx = { -- 定义域 21 { { -1, 2 }, { -1, 2 } }, 22 { { 2, 4 }, { -1, 1 } }, 23 { { 0, 4 }, { 0, 4 } }, 24 } 25 26 domv = { -- 速度限制 27 { {-0.02, 0.02}, {-0.02, 0.02} }, 28 { {-0.01, 0.01}, {-0.01, 0.01} }, 29 { {-0.02, 0.02}, {-0.02, 0.02} }, 30 } 31 32 maxfun = { -- 求解函数 33 function(x) return 4 - (x[1] * sin( 4 * pi * x[1]) - x[2] * sin(4 * pi * x[2] + pi + 1)) end , 34 function(x) return exp(x[1] - x[2] * x[2]) / (x[1] * x[1] - x[2]) + (x[1] - 3) ^ 2 end , 35 function(x) return 2 + (x[1] * x[1] + x[2] * x[2]) / 10 - (cos(2 * pi * x[1]) + cos(2 * pi * x[2])) end , 36 } 37 funidx = 0 38 -- ==================================================================== 39 40 -- ================================定义================================ 41 psonum = 80 -- 粒子种群大小 42 particles = {} -- 粒子群 43 c1 = 2 -- 学习因子(传统知识) 44 c2 = 2 -- 学习因子(社会知识) 45 domw = {0.9, 0.4} -- 惯性权重(可变) 46 w = 0.9 -- 惯性表征 47 r = 1 -- 约束因子 48 gbest = 0 -- 全局最优 49 itermax = 100 -- 迭代次数 50 -- ==================================================================== 51 52 53 -- ==============================工具函数============================== 54 -- 评价函数 55 function getfitness(y) return y end 56 57 -- 粒子群初始化 58 function psoinit() 59 -- 初始化粒子群,及每个粒子的历史最优 60 for i = 1, psonum do 61 particles[i] = { v = {}, pbest = {} } 62 for j = 1, dim do 63 particles[i][j] = domx[funidx][j][1] + (domx[funidx][j][2] - domx[funidx][j][1]) * random() 64 particles[i].v[j] = domv[funidx][j][1] + (domv[funidx][j][2] - domv[funidx][j][1]) * random() 65 end 66 particles[i].y = maxfun[funidx](particles[i]) 67 tmpfitness = getfitness(particles[i].y) 68 69 for j = 1, dim do 70 particles[i].pbest[j] = particles[i][j] 71 end 72 particles[i].ybest = particles[i].y 73 particles[i].fbest = tmpfitness 74 end 75 76 -- 更新全局最优 77 gbest = 1 78 for i = 2, psonum do 79 if particles[i].fbest > particles[gbest].fbest then gbest = i end 80 end 81 82 -- 设置惯性因子初值 83 w = domw[1] 84 end 85 86 -- 移动粒子 87 function psomove(idx) 88 for i = 1, dim do 89 -- 计算新的速度 90 particles[idx].v[i] = w * particles[idx].v[i] + 91 c1 * random() * (particles[idx].pbest[i] - particles[idx][i]) + 92 c2 * random() * (particles[gbest].pbest[i] - particles[idx][i]) 93 if particles[idx].v[i] < domv[funidx][i][1] then 94 particles[idx].v[i] = domv[funidx][i][1] 95 elseif particles[idx].v[i] > domv[funidx][i][2] then 96 particles[idx].v[i] = domv[funidx][i][2] 97 end 98 99 -- 移动粒子 100 particles[idx][i] = particles[idx][i] + r * particles[idx].v[i] 101 if particles[idx][i] < domx[funidx][i][1] then 102 particles[idx][i] = domx[funidx][i][1] 103 elseif particles[idx][i] > domx[funidx][i][2] then 104 particles[idx][i] = domx[funidx][i][2] 105 end 106 end 107 108 109 particles[idx].y = maxfun[funidx](particles[idx]) 110 tmpfitness = getfitness(particles[idx].y) 111 112 -- 更新历史最优 113 if tmpfitness > particles[idx].fbest then 114 for i = 1, dim do particles[idx].pbest[i] = particles[idx][i] end 115 particles[idx].ybest = particles[idx].y 116 particles[idx].fbest = tmpfitness 117 -- 更新全局最优 118 if tmpfitness > particles[gbest].fbest then gbest = idx end 119 end 120 end 121 122 -- ==================================================================== 123 124 125 -- ===============================主函数=============================== 126 function main(idx) 127 -- 功能选择 128 funidx = idx 129 -- 系统初始化 130 psoinit() 131 -- 惯性因子 132 deltaw = (domw[2] - domw[1]) / itermax 133 134 -- 开始迭代 135 for iter = 1, itermax do 136 -- 计算当前惯性因子 137 w = w + deltaw 138 -- 更新所有粒子位置 139 for i = 1, psonum do psomove(i) end 140 end 141 end 142 143 -- ===============================主函数=============================== 144 145 t1 = os.clock() 146 147 main(1) 148 print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2])) 149 main(2) 150 print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2])) 151 main(3) 152 print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2])) 153 154 t2 = os.clock() 155 156 print("times: ", 1000 * (t2 - t1))