python带有约束处理粒子群算法

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

算法简介

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

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

算法原理

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

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

速度$V_{i}=(v_{1},v_{2},\cdots,v_{D})$

位置$X_{i}=(x_{1},x_{2},\cdots,x_{D})$

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

$V_{i}=\omega\times V_{i}+c_{1}\times rand()\times(pbest_{i}−X_{i})+c_{2}\times rand()\times(gbest_{i}−X_{i})$

$X_{i} = X_{i} + V_{i}$

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

第一部分$V_{i}$和$X_{i}$称为惯性部分，表示上次迭代中粒子的速度、位置的影响；

第二部分$c_{1}\times rand()\times(pbest_{i}−X_{i})$称为自身认知部分，可理解为粒子当前位置与自身历史最优位置之间的距离和方向，表示粒子的动作来源于自己经验的部分；

第三部分$c_{2}\times rand()\times(gbest_{i}−X_{i})$称为群体认知部分，可理解为粒子当前位置与群体历史最优位置之间的距离和方向，反映了粒子间的协同合作和知识共享。

算法流程

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

• 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()

结果展示