粒子群算法 PSO

中规中矩, 没有ABC快

-- 粒子群算法 PSO(全局,变权重)

-- lua 相关简易操作
sin = math.sin
cos = math.cos
sqrt = math.sqrt
pi = math.pi
random = math.random
exp = math.exp
int = math.floor


-- 获得不同的随机序列
math.randomseed(os.time())


-- ==============================求解函数==============================
dim = 2                                 -- 解的维数

domx = {                                -- 定义域
    { { -1, 2 }, { -1, 2 } },
    { {  2, 4 }, { -1, 1 } },
    { {  0, 4 }, {  0, 4 } },
}

domv = {                                -- 速度限制
    { {-0.02, 0.02}, {-0.02, 0.02} },
    { {-0.01, 0.01}, {-0.01, 0.01} },
    { {-0.02, 0.02}, {-0.02, 0.02} },
}

maxfun = {                              -- 求解函数
    function(x) return 4 - (x[1] * sin( 4 * pi * x[1]) - x[2] * sin(4 * pi * x[2] + pi + 1)) end            ,
    function(x) return exp(x[1] - x[2] * x[2]) / (x[1] * x[1] - x[2]) + (x[1] - 3) ^ 2 end                  ,
    function(x) return 2 + (x[1] * x[1] + x[2] * x[2]) / 10 - (cos(2 * pi * x[1]) + cos(2 * pi * x[2])) end ,
}
funidx = 0
-- ====================================================================

-- ================================定义================================
psonum = 80                             -- 粒子种群大小
particles = {}                          -- 粒子群
c1 = 2                                  -- 学习因子(传统知识)
c2 = 2                                  -- 学习因子(社会知识)
domw = {0.9, 0.4}                       -- 惯性权重(可变)
w = 0.9                                 -- 惯性表征
r = 1                                   -- 约束因子
gbest = 0                               -- 全局最优
itermax = 100                           -- 迭代次数
-- ====================================================================


-- ==============================工具函数==============================
-- 评价函数
function getfitness(y) return y end

-- 粒子群初始化
function psoinit()
    -- 初始化粒子群,及每个粒子的历史最优
    for i = 1, psonum do
        particles[i] = { v = {}, pbest = {} }
        for j = 1, dim do
            particles[i][j] = domx[funidx][j][1] + (domx[funidx][j][2] - domx[funidx][j][1]) * random()
            particles[i].v[j] = domv[funidx][j][1] + (domv[funidx][j][2] - domv[funidx][j][1]) * random()
        end
        particles[i].y = maxfun[funidx](particles[i])
        tmpfitness = getfitness(particles[i].y)

        for j = 1, dim do
            particles[i].pbest[j] = particles[i][j]
        end
        particles[i].ybest = particles[i].y
        particles[i].fbest = tmpfitness
    end

    -- 更新全局最优
    gbest = 1
    for i = 2, psonum do
        if particles[i].fbest > particles[gbest].fbest then gbest = i end
    end

    -- 设置惯性因子初值
    w = domw[1]
end

-- 移动粒子
function psomove(idx)
    for i = 1, dim do
        -- 计算新的速度
        particles[idx].v[i] = w * particles[idx].v[i] + 
            c1 * random() * (particles[idx].pbest[i] - particles[idx][i]) +
            c2 * random() * (particles[gbest].pbest[i] - particles[idx][i])
        if particles[idx].v[i] < domv[funidx][i][1] then
            particles[idx].v[i] = domv[funidx][i][1]
        elseif particles[idx].v[i] > domv[funidx][i][2] then
            particles[idx].v[i] = domv[funidx][i][2]
        end
        
        -- 移动粒子
        particles[idx][i] = particles[idx][i] + r * particles[idx].v[i]
        if particles[idx][i] < domx[funidx][i][1] then
            particles[idx][i] = domx[funidx][i][1]
        elseif particles[idx][i] > domx[funidx][i][2] then
            particles[idx][i] = domx[funidx][i][2]
        end
    end
    
    
    particles[idx].y = maxfun[funidx](particles[idx])
    tmpfitness = getfitness(particles[idx].y)

    -- 更新历史最优
    if tmpfitness > particles[idx].fbest then
        for i = 1, dim do particles[idx].pbest[i] = particles[idx][i] end
        particles[idx].ybest = particles[idx].y
        particles[idx].fbest = tmpfitness
        -- 更新全局最优
        if tmpfitness > particles[gbest].fbest then gbest = idx end
    end
end

-- ====================================================================


-- ===============================主函数===============================
function main(idx)
    -- 功能选择
    funidx = idx
    -- 系统初始化
    psoinit()
    -- 惯性因子
    deltaw = (domw[2] - domw[1]) / itermax
    
    -- 开始迭代
    for iter = 1, itermax do
        -- 计算当前惯性因子
        w = w + deltaw
        -- 更新所有粒子位置
        for i = 1, psonum do psomove(i) end
    end
end

-- ===============================主函数===============================

t1 = os.clock()

main(1)
print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2]))
main(2)
print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2]))
main(3)
print(string.format("函数值为: %.8f \t解为: (%.3f, %.3f)", particles[gbest].y, particles[gbest].pbest[1], particles[gbest].pbest[2]))

t2 = os.clock()

print("times: ", 1000 * (t2 - t1))