優(yōu)化算法庫DEAP的粒子群優(yōu)化算法（PSO）示例代碼分析

之前使用優(yōu)化算法庫DEAP做遺傳算法(Genetic Optimization)，最近在研究粒子群優(yōu)化算法（PSO）發(fā)現(xiàn)DEAP也提供支持；對原代碼做了分析，之后會嘗試用于交易策略參數(shù)優(yōu)化。

首先關(guān)于 粒子群優(yōu)化算法（PSO） ，這個算法的介紹很多，我看了不少，感覺下文介紹的是比較清楚的，可以參考。

https://www.cnblogs.com/21207-iHome/p/6062535.html

其實原理很簡單，看看這個圖就可以，可以理解為一群蜜蜂去公園采蜜，每一輪蜜蜂都交流，看那個蜜蜂聞到味道最濃，大家都往那邊去，然后每輪更新，最后聚集到最濃的地方，就是最優(yōu)位置。

這樣，PSO把優(yōu)化的對象設置為位置點，比如兩個參數(shù)要優(yōu)化，就是（x1,x2)； 如果有N個參數(shù)要優(yōu)化，那么就是N維的位置點，可以用list放置；和遺傳算法GA差別就是，引入了速度v，可以理解參數(shù)的變動的速率；同時還有一個集群最優(yōu)參數(shù)，就是每一輪所有粒子共享信息后，要去聚集位置，每輪更新。

其他引用鏈接講的很詳細。主要過程就是如下

• 粒子群算法流程 

第1步   在初始化范圍內(nèi)，對粒子群進行隨機初始化，包括隨機位置和速度

第2步   計算每個粒子的適應值

第3步   更新粒子個體的歷史最優(yōu)位置

第4步   更新粒子群體的歷史最優(yōu)位置

第5步   更新粒子的速度和位置

第6步   若未達到終止條件，則轉(zhuǎn)第2步

　　粒子群算法流程圖如下：pbest是個體的最優(yōu)質(zhì)，gbest是粒子群中最優(yōu)的粒子

DEAP的示例代碼，其實也是按照這個邏輯。

1，首先，通過類定義方法creator.create，定義兩個類， 可以看我之前將GA示例代碼介紹。

第一個類FitnessMax，是用于定義按最大值優(yōu)化；如果weights = (-1.0,) 就是選取最小值。

第二個類Particle，是的定義粒子；可以看到第一個list是這個 Particle 繼承于list，list用來存放空間位置，同時有一系列參數(shù)。

其中FitnessMax定義優(yōu)化方向，speed也是一個list，用來存放每個參數(shù)變化步進，smin和smax是速度的上下限，best是歷史最優(yōu)參數(shù)。

creator.create("FitnessMax", base.Fitness, weights=(1.0,))
creator.create("Particle", list, fitness=creator.FitnessMax, speed=list, 
    smin=None, smax=None, best=None)

2，然后，定義粒子初始化方法generate，其中size是粒子的維度，比如有10個參數(shù)要優(yōu)化，那么size就是10；其中pmin，pmax是創(chuàng)建時候可能的位置上下限；smin，smax是創(chuàng)建時候的速度上線。這個方法就生成了一個粒子，包括初始時候位置，速度，和速度上下限；只有best最優(yōu)參數(shù)還沒有計算。

def generate(size, pmin, pmax, smin, smax):
    part = creator.Particle(random.uniform(pmin, pmax) for _ in range(size)) 
    part.speed = [random.uniform(smin, smax) for _ in range(size)]
    part.smin = smin
    part.smax = smax
    return part

3，然后，在定義粒子更新方法updateParticle，這里代碼有就是為了實現(xiàn)下面兩個公式，更新速度，并且根據(jù)速度和位置求出更新位置。

輸入best參數(shù)是粒子群所有粒子中的最優(yōu)，part是更新的粒子，ph2，ph3就是c1r1隨機函數(shù)和加速度常數(shù)積。唯一區(qū)別代碼沒有使用慣性權(quán)重。

def updateParticle(part, best, phi1, phi2):
        #求出新速度，如速度更新公式
    u1 = (random.uniform(0, phi1) for _ in range(len(part)))
    u2 = (random.uniform(0, phi2) for _ in range(len(part)))
    v_u1 = map(operator.mul, u1, map(operator.sub, part.best, part))
    v_u2 = map(operator.mul, u2, map(operator.sub, best, part))
    part.speed = list(map(operator.add, part.speed, map(operator.add, v_u1, v_u2)))
    #求出新位置，這里如速度超過smin，smax，則使用smin
    for i, speed in enumerate(part.speed):
        if speed < part.smin:
            part.speed[i] = part.smin
        elif speed > part.smax:
            part.speed[i] = part.smax
    part[:] = list(map(operator.add, part, part.speed))

4、然后就把上面這些定義綁定在toolbox里面，這樣就方便批量調(diào)用。

    - particle：綁定粒子生成方法，粒子位置有兩位，就像[x,y]，位置和速度最大最小值

    -  population：生成多個粒子

    - update：粒子更新方法

    -   evaluate： 計算最優(yōu)值，這里使用了benchmarks.h2，我把源代碼列如下。要求返回值最大，可以看到返回值是num / denum， 只要 denum最小就可以，那么individual[0]是8.6998， individual[1]是6.7665，就是最小了。

toolbox = base.Toolbox()
toolbox.register("particle", generate, size=2, pmin=-6, pmax=6, smin=-3, smax=3)
toolbox.register("population", tools.initRepeat, list, toolbox.particle)
toolbox.register("update", updateParticle, phi1=2.0, phi2=2.0)
toolbox.register("evaluate", benchmarks.h2)
def h2(individual):
    num = (sin(individual[0] - individual[1] / 8))**2 + (sin(individual[1] + individual[0] / 8))**2
    denum = ((individual[0] - 8.6998)**2 + (individual[1] - 6.7665)**2)**0.5 + 1
    return num / denum,

5、后面就是執(zhí)行了，我把源代碼的log代碼刪除了，比較簡單，我直接注釋了。

def main():
    pop = toolbox.population(n=5) #粒子群有5個粒子
    GEN = 1000 #更新一千次
    best = None 
    for g in range(GEN):
        for part in pop: #每次更新，計算粒子群中最優(yōu)參數(shù)，并把最優(yōu)值寫入best
            part.fitness.values = toolbox.evaluate(part) 
            if not part.best or part.best.fitness < part.fitness:
                part.best = creator.Particle(part)
                part.best.fitness.values = part.fitness.values
            if not best or best.fitness < part.fitness:
                best = creator.Particle(part)
                best.fitness.values = part.fitness.values
        for part in pop: #更新粒子位置
            toolbox.update(part, best)
        # Gather all the fitnesses in one list and print the stats
    return pop, best

然后返回運行，得到的最優(yōu)位置如下，和準確的差不多。

[8.720808209484728, 6.765281780793948]

源代碼鏈接如下

https://github.com/DEAP/deap/blob/82f774d9be6bad4b9d88272ba70ed6f1fca39fcf/examples/pso/basic.py