优化算法|粒子群优化算法及python代码实现

作者简介:本人擅长运筹优化建模及算法设计,包括各类车辆路径问题、生产车间调度、二三维装箱问题,熟悉CPLEX和gurobi求解器
微信公众号:运筹优化与学习

若有运筹优化建模及算法定制需求,欢迎联系我们私聊沟通

算法简介

粒子群算法(Particle Swarm Optimization, PSO)的思想源于对鸟群觅食行为的研究,其核心思想是通过群体中个体之间的协作和信息共享来寻找最优解。

相较于遗传算法,粒子群算法具有收敛速度快、参数少、算法简单易实现的优点(对高维度优化问题,比遗传算法更快收敛于最优解),但是同样存在陷入局部最优解的问题。

算法原理

粒子群算法通过设计粒子来模拟鸟群中的鸟,粒子具有两个属性:速度和位置,其中速度表示粒子下一步迭代时移动的方向和距离,位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解,并将其记为当前个体极值,并将个体极值与整个粒子群里的其他粒子共享,找到最优的那个个体极值作为整个粒子群的当前全局最优解,粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。

假设在 D D D维搜索空间 S S S中,有 N N N个粒子,每个粒子代表一个解,则第 i i i个粒子的速度、位置表达式如下:

速度 V i = ( v 1 , v 2 , ⋯   , v D ) V_{i}=(v_{1},v_{2},\cdots,v_{D}) Vi=(v1,v2,,vD)

位置 X i = ( x 1 , x 2 , ⋯   , x D ) X_{i}=(x_{1},x_{2},\cdots,x_{D}) Xi=(x1,x2,,xD)

在每一次的迭代中,粒子通过跟踪两个“极值”(pbest,gbest)来更新自己的速度和位置,更新表达式如下:

V i = ω × V i + c 1 × r a n d ( ) × ( p b e s t i − X i ) + c 2 × r a n d ( ) × ( g b e s t i − X i ) V_{i}=\omega\times V_{i}+c_{1}\times rand()\times(pbest_{i}−X_{i})+c_{2}\times rand()\times(gbest_{i}−X_{i}) Vi=ω×Vi+c1×rand()×(pbestiXi)+c2×rand()×(gbestiXi)

X i = X i + V i X_{i} = X_{i} + V_{i} Xi=Xi+Vi

式中, ω \omega ω表示惯性因子,现有研究表明采用动态惯性因子能够获得更好的寻优结果,目前最常用的为线性递减权值(Linear Decreasing Weight,LDW)策略。 c 1 c_1 c1为个体学习因子, c 2 c_2 c2为群体学习因子。pbest表示当前粒子搜索最好解,gbest表示群体中搜索最好解。

第一部分 V i V_{i} Vi X i X_{i} Xi称为惯性部分,表示上次迭代中粒子的速度、位置的影响;

第二部分 c 1 × r a n d ( ) × ( p b e s t i − X i ) c_{1}\times rand()\times(pbest_{i}−X_{i}) c1×rand()×(pbestiXi)称为自身认知部分,可理解为粒子当前位置与自身历史最优位置之间的距离和方向,表示粒子的动作来源于自己经验的部分;

第三部分 c 2 × r a n d ( ) × ( g b e s t i − X i ) c_{2}\times rand()\times(gbest_{i}−X_{i}) c2×rand()×(gbestiXi)称为群体认知部分,可理解为粒子当前位置与群体历史最优位置之间的距离和方向,反映了粒子间的协同合作和知识共享。

算法流程

粒子群优化算法的步骤如下:

  • Step 1:初始化粒子速度、位置,产生初始种群
  • Step 2:计算当前种群中每个解的适应度(目标函数值)
  • Step 3:更新个体最优解与群体最优解
  • Step 4:利用速度、位置更新公式,更新当前解集
  • Step 5:计算当前种群里每个解的适应度(目标函数值)
  • Step 6:更新个体最优解与群体最优解,若满足迭代终止条件时,跳出迭代环节,否则重复执行Step 4-6
  • Step 7:输出搜索到的最优解

单目标代码展示

定义粒子类

class particle_single():
    '''
    position -> 粒子位置
    velocity -> 粒子速度
    best_p -> 全局最优解
    obj_function -> 目标函数
    constraints -> 约束
    obj_value -> 目标函数值
    '''
    
    def __init__(self,obj_func,attribute_number,constr=[],vmax=np.array(np.nan),l_bound=np.array(np.nan),u_bound=np.array(np.nan),integer=np.array(np.nan),position=np.array(np.nan),velocity=np.array(np.nan)):
        self.obj_function = obj_func
        if type(constr)!=list:
            constr = [constr]
        self.constraints = constr
        self.att = attribute_number
        if np.all(np.isnan(position))==False:
            self.position = position
        else:
            try:
                if attribute_number!=1:
                    init_pos = []
                    for i in range(self.att):
                        if integer[i]==False:
                            init_pos.append(random.uniform(l_bound[i],u_bound[i]))
                        else:
                            init_pos.append(random.randint(l_bound[i],u_bound[i]))
                    self.position = np.array(init_pos)
                else:
                    if integer==False:
                        self.position= random.uniform(l_bound,u_bound)
                    else:
                        self.position= random.randint(l_bound,u_bound)
            except:
                print('We need lower and upper bound for init position')
        
        self.obj_value = self.calc_obj_value()
        if np.all(np.isnan(velocity))==False:
            self.velocity = velocity
        else:
            try:
                if attribute_number!=1:
                    self.velocity = np.array([random.uniform(-vmax[i],vmax[i]) for i in range(self.att)])
                else:
                    self.velocity = random.uniform(-vmax,vmax)
            except:
                print('we need an vmax for init velocity')
        self.best_p=np.nan

多目标代码展示

测试函数

定义粒子类

class particle_multi(particle_single):
    '''
    position -> 粒子位置
    velocity -> 粒子速度
    best_p -> best particle in the past gets just set by first_set_best
    obj_functions -> 目标函数
    constraints -> 约束条件
    obj_values -> 目标函数值
    S -> 当前解支配解集
    n -> 支配解数量
    distance -> 拥挤距离
    rank -> 支配排序
    '''
    
    def __init__(self,obj_func,attribute_number,constr=[],vmax=np.array(np.nan),l_bound=np.array(np.nan),u_bound=np.array(np.nan),integer=np.array(np.nan),position=np.array(np.nan),velocity=np.array(np.nan)):
        self.obj_functions = obj_func
        self.constraints=constr
        self.att = attribute_number
        if np.all(np.isnan(position))==False:
            self.position=position
        else:
            try:
                if attribute_number!=1:
                    init_pos = []
                    for i in range(self.att):
                        if integer[i]==False:
                            init_pos.append(random.uniform(l_bound[i],u_bound[i]))
                        else:
                            init_pos.append(random.randint(l_bound[i],u_bound[i]))
                    self.position = np.array(init_pos)
                else:
                    if integer==False:
                        self.position= random.uniform(l_bound,u_bound)
                    else:
                        self.position= random.randint(l_bound,u_bound)
            except:
                print('We need lower and upper bound for init position')
        
        self.obj_values =self.calc_obj_value()
        
        if np.all(np.isnan(velocity))==False:
            self.velocity = velocity
        else:
            try:
                if attribute_number!=1:
                    self.velocity = np.array([random.uniform(-vmax[i],vmax[i]) for i in range(self.att)])
                else:
                    self.velocity = random.uniform(-vmax,vmax)
            except:
                print('we need an vmax for init velocity')
        self.best_p=np.nan
        self.S = []
        self.n = np.nan
        self.rank = np.nan
        self.distance = np.nan

算法主循环

def runner(self,Iterations, time_termination):
        t0 = time.time()
        for i in range(Iterations):  # 迭代次数上限
            if time_termination != -1 and time.time()-t0 > time_termination: # 求解时间上限
                break
            if self.multi:
                self.comp_swarm=[]
                self.comp_swarm.append(self.g_best)
            for part in self.swarm:
                # 更新速度
                r1 = random.random()
                r2 = random.random()
                new_v = self.v_weight*part.velocity + self.c_param*r1*(part.best_p.position-part.position) + self.s_param*r2*(self.g_best.position-part.position)
                # control vmax
                new_v = np.array([new_v[i] if new_v[i]>-self.vmax[i] else -self.vmax[i] for i in range(len(new_v))])
                new_v = np.array([new_v[i] if new_v[i]<self.vmax[i] else self.vmax[i] for i in range(len(new_v))])                
                
                # 更新位置
                new_p = part.position + new_v
                for i in range(len(new_p)):
                    if self.integer[i]:
                        new_p[i] = int(new_p[i])
                
                # 检验边界条件
                new_p = np.array([new_p[i] if new_p[i]>self.l_bound[i] else self.u_bound[i] - abs(self.l_bound[i]-new_p[i])%(self.u_bound[i]-self.l_bound[i]) for i in range(len(new_p))])
                new_p = np.array([new_p[i] if new_p[i]<self.u_bound[i] else self.l_bound[i] + abs(self.u_bound[i]-new_p[i])%(self.u_bound[i]-self.l_bound[i]) for i in range(len(new_p))])               
                
                part.set_velocity(new_v)
                part.set_position(new_p)
                
                if self.multi:
                    self.comp_swarm.append(part)
                    self.comp_swarm.append(part.best_p)
                
            if self.multi:
                self.non_dom_sort()  # 非支配排序
                self.g_best = copy.deepcopy(self.comp_swarm[0])
                j=1
                new_swarm = self.comp_swarm[1:-1:2]
                self.swarm = copy.deepcopy(new_swarm)
                j+=1
                
            for part in self.swarm:
                if self.multi:
                    part.best_p = copy.deepcopy(self.comp_swarm[j])
                    j+=2
                if part.compare(self.g_best):
                    self.g_best = copy.deepcopy(part)
                # 更新当前最优解
                part.compare_p_best()

结果展示

若有运筹优化建模及算法定制需求,欢迎联系我们私聊沟通

  • 22
    点赞
  • 40
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
以下是Python实现粒子优化算法的示例代码: ```python import random import math class Particle: def __init__(self, x, y): self.x = x self.y = y self.vx = 0 self.vy = 0 self.bestx = x self.besty = y self.bestscore = float('inf') self.score = float('inf') def update_velocity(self, bestx, besty, w, c1, c2): self.vx = w * self.vx + c1 * random.random() * (self.bestx - self.x) + c2 * random.random() * (bestx - self.x) self.vy = w * self.vy + c1 * random.random() * (self.besty - self.y) + c2 * random.random() * (besty - self.y) def update_position(self): self.x += self.vx self.y += self.vy def evaluate(self, func): self.score = func(self.x, self.y) if self.score < self.bestscore: self.bestscore = self.score self.bestx = self.x self.besty = self.y def particle_swarm(func, bounds, n_particles, n_iter, w, c1, c2): particles = [Particle(random.uniform(bounds[0][0], bounds[1][0]), random.uniform(bounds[0][1], bounds[1][1])) for i in range(n_particles)] bestscore = float('inf') bestx = None besty = None for i in range(n_iter): for particle in particles: particle.evaluate(func) if particle.score < bestscore: bestscore = particle.score bestx = particle.x besty = particle.y for particle in particles: particle.update_velocity(bestx, besty, w, c1, c2) particle.update_position() particle.x = max(min(particle.x, bounds[1][0]), bounds[0][0]) particle.y = max(min(particle.y, bounds[1][1]), bounds[0][1]) return (bestx, besty, bestscore) ``` 该代码定义了一个 `Particle` 类来表示每个粒子,包含了位置、速度、历史最优位置和历史最优得分等属性,以及更新速度、更新位置和评估得分等方法。 `particle_swarm` 函数实现粒子优化算法,接受函数、变量边界、粒子数、迭代次数、惯性权重、个体学习因子和群体学习因子等参数。在函数内部,先初始化了一组粒子,然后在每次迭代中,先评估每个粒子的得分,更新全局最优解,然后更新每个粒子的速度和位置,并确保它们在变量边界范围内。最后返回全局最优解的位置和得分。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

eternal1995

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值