之前使用优化算法库DEAP做遗传算法(Genetic Optimization),最近在研究粒子群优化算法(PSO)发现DEAP也提供支持;对原代码做了分析,之后会尝试用于交易策略参数优化。
首先关于
粒子群优化算法(PSO)
,这个算法的介绍很多,我看了不少,感觉下文介绍的是比较清楚的,可以参考。
其实原理很简单,看看这个图就可以,可以理解为一群蜜蜂去公园采蜜,每一轮蜜蜂都交流,看那个蜜蜂闻到味道最浓,大家都往那边去,然后每轮更新,最后聚集到最浓的地方,就是最优位置。
这样,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.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.best = creator.Particle(part)
part.best.fitness.values = part.fitness.values
if not best or best.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