富联注册: 叶菜类  根茎类  果实类  主食厨房 
全国服务热线:0898-08980898
富联资讯
当前位置: 首页 > 富联资讯
粒子群算法求解带约束优化问题
添加时间:2024-06-04
  

我们在生活中经常课件带约束的优化问题,比如蛋糕的分配,店铺的选址,如何以最低费用制定旅行计划等。这些优化问题,都可以表示如下形式:

\\min_x f(x), \\ x \\in R \\\\ g_i(x) \\leq 0, \\ i=1, .. m

其中, m是约束的个数,f(x)是目标函数, g_i(x)表示第i个约束方程。大于0的约束可以表示成

-g_i(x) \\leq 0

等于0的约束可以写为:

-g_i(x) \\leq 0 \\\\ g_i(x) \\leq 0

通常地,我们针对该问题可以通过将数学公式编写为代码,求解器进行求解。如果遇到问题规模较大,决策变量和约束较多,启发式搜索的方法可以被采用。本文介绍用粒子群算法求解带约束问题的方法。

粒子群算法,英文为Particle Swarm Optimization,简称为PSO,又称为微粒群算法。是由 J. Kennedy 和 R. C. Eberhart 等于1995年开发的一种演化计算技术,来源于对一个简化社会模型的模拟。

PSO 算法最初是为了图形化地模拟鸟群优美而不可预测的运动。而通过对动物社会行为的观察,发现在群体中对信息的社会共享提供一个演化的优势,并以此作为开发算法的基础。通过加入近邻的速度匹配、并考虑了多维搜索和根据距离的加速,形成了 PSO 的最初版本。之后引入了惯性权重w来更好的控制开发(exploitation)和探索(exploration),形成了标准版本。为了提高粒群算法的性能和实用性,中山大学、(英国)格拉斯哥大学等又开发了自适应(Adaptive PSO)版本和离散(discrete)版本。[2]

粒子群算法的更新公式如下:

\\begin{align}V_i^{k+1}&:=w V_i^k + c_1 r_{i1}^k (P_i^k - X_i^k) + c_2 r_{i2}^k (P_g^k - X_i^k)  \\\\ X_i^{k+1}&:=X_i^{k}+ V_i^{k+1}\\end{align}

其中,

  • V_i^{k+1} : 表示第 i 个例子,第 k+1 次迭代的速度,初始化速度为0
  • w : 是惯性权重,如前所述用于控制exploitation和exploration
  • c_1 :称为认知常数,实现自我认知。通常取值为2
  • c_2 :称为社会常数,该部分模拟与社会种群的交互,通常取值也为2
  • r_{i1}^kr_{i2}^k :两个正态分布取值的随机数,范围为[0,1]
  • P_i^k :是截止到第 k 次迭代,第i个粒子获得最优目标/fitness值对应的解
  • P_g^k :是截止到第 k 次迭代,所有粒子中获得最优目标/fitness值对应的解
  • X_i^k :是第 k 次迭代,所有粒子对应的解,矩阵维度为(粒子数量,解维度)。初始解可以初始化为在所有解维度中允许取值范围的并集。比如 0\\leq x_1 \\leq105 \\leq x_2 \\leq 12 。初始化可以为0-12的正态分布,这样可以加速搜索。

公式1对应粒子速度的迭代,公式2更新当前解。

上述介绍了粒子群的主要部分。作为演化算法中的一员,粒子群算法的执行还需要有一个重要的部分就是fitness函数,也可以成为解的评价函数。fitness函数可以为目标值 f(x) ,也可以为其它的函数,它的作用主要就是来给解进行评分。

整体的粒子群算法求解流程如下:

  1. 初始化参数,粒子群,速度等,保存当前的P_g P_i 和最优目标值
  2. 计算当前粒子群各粒子的fitness分数,更新P_g P_i 和最优目标值
  3. 更新速度和解
  4. 判断是否达到最大迭代数,如果达到则退出,如无则返回第2步

根据前面的粒子群算法介绍,粒子群算法中没有涉及求解带约束的优化问题。那么约束如何加入呢?下面会介绍粒子群求解带约束优化问题的理论基础和实践。

理论基础

有些同学马上能够反映是用拉格朗日乘子法给目标函数加入约束的惩罚项。思路是类似的,我们把目标函数变为:

\\begin{align*}F(x) &=f(x) + h(k)H(x) \\\\ h(k) &=\\sqrt{k}\\ \	ext{or}\\ k\\sqrt{k}\\\\ H(x) &=\\sum_{i=1}^m \	heta (q_i(x)) q_i(x_i)^{\\gamma(q_i(x_i))}\\end{align*}

  • f(x) :是原来的目标函数
  • h(x) :动态更新的惩罚系数,与迭代次数相关。通常可使用 k\\sqrt{k} 。如果问题规模更小,可以使用 \\sqrt{k}
  • H(x) :是约束惩罚项

其中,

  • q_i(x) 是相对约束惩罚函数
  • \	heta(q_i(x)) 是分段赋值函数
  • \\gamma(q_i(x)) 是惩罚指数

这个三个值是依据具体问题来确定的,以下提供一个相对通用的参考。

q_i(x)=\\max(0, g_i(x))\\\\ \	heta(q_i(x))=\\begin{cases}10 & q_i(x) < 0.001 \\\\ 20 & 0.001 \\leq q_i(x) \\leq 0.1 \\\\ 100 & 0.1 < q_i(x) \\leq 1 \\\\  300 & o.w.  \\end{cases}\\\\ \\gamma(q_i(x))=\\begin{cases}1 & q_i(x) < 1 \\\\  2 & o.w.  \\end{cases}

实践

在具体实践中,为了防止速度过大而导致更新的解爆炸(explosion),我们会在迭代公式中加入速度限制因子 0<\\chi<1 , 比如 \\chi 可以取值为0.7.

\\begin{align}V_i^{k+1}&:=\\chi (w V_i^k + c_1 r_{i1}^k (P_i^k - X_i^k) + c_2 r_{i2}^k (P_g^k - X_i^k)) \\\\ X_i^{k+1}&:=X_i^{k}+ V_i^{k+1}\\end{align}

同时,我们设置最大速度 V_{max} ,使得速度控制在 [-V_{max}, V_{max}] 范围内。

该步骤可以类比深度学习中为了防止梯度爆炸,做梯度裁剪(gradient clip)。

在迭代过程中,类似于SGD中学习率逐渐递减,粒子群算法中的惯性权重也是随着迭代次数的增加逐渐减少。这里惯性权重初始化为1.2,逐渐减少到0.1。

依照上面的步骤进行代码编写成solver,详见如下链接:

Github 链接:

jingw2/solver

使用示例:

假如我们要求解这个优化问题

\\begin{align*}\\min f(x) &=(x_1 - 2)^2 + (x_2 - 1)^2 \\\\ s.t. \\ x_1 &=2x_2 - 1 \\\\ x_1^2 /4 &+ x_2^2 - 1 \\leq 0 \\\\ \\end{align*}

变为标准格式 g_i(x) \\leq 0

\\begin{align*}\\min f(x) &=(x_1 - 2)^2 + (x_2 - 1)^2 \\\\ s.t. \\ x_1 &- 2x_2 + 1 \\leq 0 \\\\ - x_1 &+ 2x_2 - 1 \\leq 0 \\\\ x_1^2 /4 &+ x_2^2 - 1 \\leq 0 \\\\ \\end{align*}

代码实现:

import psoco
import math 

def objective(x):
    '''create objectives based on inputs x as 2D array'''
    return (x[:, 0] - 2) ** 2 + (x[:, 1] - 1) ** 2 

def constraints1(x):
    '''create constraint1 based on inputs x as 2D array'''
    return x[:, 0] - 2 * x[:, 1] + 1 

def constraints2(x):
    '''create constraint2 based on inputs x as 2D array'''
    return - (x[:, 0] - 2 * x[:, 1] + 1)

def constraints3(x):
    '''create constraint3 based on inputs x as 2D array'''
    return x[:, 0] ** 2 / 4. + x[:, 1] ** 2 - 1

def new_penalty_func(k):
    '''Easy Problem can use \\sqrt{k}'''
    return math.sqrt(k)
    
constraints = [constraints1, constraints2, constraints3]
num_runs = 10
# random parameters lead to variations, so run several time to get mean
for _ in range(num_runs):
    pso = psoco.PSOCO(sol_size=2, fitness=objective, constraints=constraints)
    pso.h = new_penalty_func
    pso.init_Population(low=0, high=1) # x并集的上下限,默认为0和1
    pso.solve()
    # best solutions
    x = pso.gbest.reshape((1, -1))
    # best fit which is not objective
    fit = pso.fit

可以得到最优的目标值1.39左右。

目前该package已经上传Pypi,详细链接请点击如下链接。可用pip进行安装。

psoco

[1]Particle Swarm Optimization Method for Constrained Optimization Problems: cs.cinvestav.mx/~constr

[2]维基百科:zh.wikipedia.org/wiki/%

[3]blog.csdn.net/google198

0898-08980898
手机:
13966668888
邮箱:
admin@youweb.com
电话:
0898-08980898
地址:
海南省海口市

平台注册入口