優化演算法庫DEAP的粒子群優化演算法（PSO）示例程式碼分析

之前使用優化演算法庫DEAP做遺傳演算法(Genetic Optimization)，最近在研究粒子群優化演算法（PSO）發現DEAP也提供支援；對原始碼做了分析，之後會嘗試用於交易策略引數優化。

首先關於 粒子群優化演算法（PSO） ，這個演算法的介紹很多，我看了不少，感覺下文介紹的是比較清楚的，可以參考。

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

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

這樣，PSO把優化的物件設定為位置點，比如兩個引數要優化，就是（x1,x2)； 如果有N個引數要優化，那麼就是N維的位置點，可以用list放置；和遺傳演算法GA差別就是，引入了速度v，可以理解引數的變動的速率；同時還有一個叢集最優引數，就是每一輪所有粒子共享資訊後，要去聚集位置，每輪更新。

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

• 粒子群演算法流程 

第1步   在初始化範圍內，對粒子群進行隨機初始化，包括隨機位置和速度

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

第3步   更新粒子個體的歷史最優位置

第4步   更新粒子群體的歷史最優位置

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

第6步   若未達到終止條件，則轉第2步

　　粒子群演算法流程圖如下：pbest是個體的最優質，gbest是粒子群中最優的粒子

 

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

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

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

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

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

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個引數要優化，那麼size就是10；其中pmin，pmax是建立時候可能的位置上下限；smin，smax是建立時候的速度上線。這個方法就生成了一個粒子，包括初始時候位置，速度，和速度上下限；只有best最優引數還沒有計算。

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，這裡程式碼有就是為了實現下面兩個公式，更新速度，並且根據速度和位置求出更新位置。

輸入best引數是粒子群所有粒子中的最優，part是更新的粒子，ph1，ph2就是c1r1隨機函式和加速度常數積。唯一區別程式碼沒有使用慣性權重。

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裡面，這樣就方便批量呼叫。

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

    -  population：生成多個粒子

    - update：粒子更新方法

    -   evaluate： 計算最優值，這裡使用了benchmarks.h1，我把原始碼列如下。要求返回值最大，可以看到返回值是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.h1)
def h1(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、後面就是執行了，我把原始碼的log程式碼刪除了，比較簡單，我直接註釋了。

def main():
    pop = toolbox.population(n=5) #粒子群有5個粒子
    GEN = 1000 #更新一千次
    best = None 
    for g in range(GEN):
        for part in pop: #每次更新，計算粒子群中最優引數，並把最優值寫入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

然後返回執行，得到的最優位置如下，和準確的差不多。

[8.720808209484728, 6.765281780793948]

原始碼連結如下

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

來自 “ ITPUB部落格 ” ，連結：http://blog.itpub.net/22259926/viewspace-2644230/，如需轉載，請註明出處，否則將追究法律責任。