Unity开发者联盟收藏本站
发帖教程新人必读切换到宽版

 找回密码
 立即注册
Unity开发者联盟»论坛 › 3D Answers 技术讨论 › 技术讨论 Technology Answers › 智能优化算法学习-粒子群算法（Particle Swarm Optimiza ...
返回列表 发新帖
查看: 141|回复: 0

智能优化算法学习-粒子群算法（Particle Swarm Optimization）-2解决约束优化问题(附python代码)

[复制链接]
发表于 2024-7-15 17:44 | 显示全部楼层 |阅读模式
想详细了解PSO的理论部门请移步到：
智能优化算法学习-粒子群算法（Particle Swarm Optimization）-1解决无约束优化问题(附python代码)
这一次的重点主要是学习了废弛约束的方式，PSO的算法只是求解的一种途径，还可以换其他的智能优化算法来求解。（下午在健身得时候，师姐问了一下lower bound technique是否属于精确性算法，然后在兴禄学长那边得到了更全面的认识，十分感激！
精确不精确，主要看最后能不能guarantee optimality，能不能主要给充沛的时间，就能在有限的步数内收敛到最优解，如果可以则属于精确性算法，只能给lower bound和upper bound，但是不能保证收敛到精确解，这种叫做近似算法，这里需要注意这个“足够的时间”，可以是几十年，甚至是几亿年）
思路是将约束作为惩罚项添加到方针函数中（类似Interior Point Method，Lagrangian Relaxation），将模型转化为无约束的模型，进而求解（例如，操作牛顿法求解无约束优化问题）。下面操作一个例题来进行说明：
假设有如下带约束的优化问题：
minf(x)      
s.t.:g_{i}(x) \leq0, i=1, 2,....,l                               ①
h_{j}(x)=0,j=0,1,....,m           ②
约束中既含有不等式约束和等式约束，当一个解x不在可行空间内（即不满足约束条件），则会对方针函数增加惩罚，这样就把有约束优化问题转化为无约束优化问题，转化后的方针函数如下：
F(x,\sigma)=f(x)+\sigma P(x)
此中 \sigma 是惩罚因子， P(x) 是整体惩罚项，对约束的措置：
1) 不等式约束：对不不等式①我们可以直接放入方针函数中与惩罚因子相乘，不等式约束可以构造为 \leq 0的形式，直接在计算方针是进行惩罚。
这里对于不等式约束 g(x)\leq0 的惩罚项 e_{i}=max(0,-g_{i}(x)) ,即当 g_{i}(x)\leq0 时（满足约束），惩罚项为0，否则等于 -g_{i}(x) 。
2) 等式约束：因为是要作为惩罚项放到方针函数中，最好在不满足约束时可以直接对方针函数造成影响。这里选用一个简单的措置方式——设定一个等式约束容忍度值 \varepsilon ，将其转化成一个不等式约束，则新的不等式约束则可以暗示为 |h_{j}(x)|\leq\varepsilon\Rightarrow|h_{j}(x)|-\varepsilon\leq0 ,然后同上措置 e_{j}(x)=max(0, |h_{j}(x)|-\varepsilon)
综上 P(x)=\sum_{i=1}^{l}{e_{i}(x)}+\sum_{j=1}^{m}{e_{j}(x)}
=\sum_{i=1}^{l}max({0,-g_{i}(x))}+\sum_{j=i+1}^{l}max({0,|h_{j}(x)|-\varepsilon)}
由于约束条件之间存在差异，某些约束条件可能对个体违反约束程度起着决定性感化，为了平衡这种差异，对惩罚项做归一化措置，为了便利起见，在推导公式中不区分等式约束与不等式约束，在实际措置时再做区分即可：
L_{j}=\sum_{i=1}^{N}{e_{j}(x_{i})}/\sum_{j=1}^{m}\sum_{i=1}^{N}{e_{j}(x_{i})}, j=1,2,...,m
此中 L_{j} 暗示每个约束条件违背程度，m为约束条件的个数。公式中分子的意思是对于每个粒子 x_{i} 计算违反第j个约束的惩罚项，分母则暗示对每个粒子 x_{i} 计算违反每个约束的惩罚项再作求和。故L_{j} 可以视为第j个惩罚项的权重。
最后得到粒子 x_{i} 的整体惩罚项 P(x) 的计算公式：
P(x_{i})=L_{j}e_{j}(x_{i})
\sum_{j=1}^{m}{L_{j}}=1
在粒子群算法中，每一步迭代城市更新Pbest和Gbest，虽然可以将有约束问题转成无约束问题进行迭代求解，但是在求解过程中依然会存在不满足约束条件的情况，因此需要必然的法则来斗劲两个粒子的优劣：
①如果两个例子 x_{i}和x_{j} 都可行，则斗劲器适应度函数值 f(x_{i}) 和 f(x_{j})  ，值小为优。
②当两个粒子 x_{i}和x_{j} 都不成行，则斗劲惩罚项 P(x_{i}) 和 P(x_{j}) ，违反约束程度小的粒子为优。
③当粒子 x_{i} 可行，而 x_{j} 不成行，则选可行解。
<hr/>同样的，我们还是选用Rastrigin函数作为测试函数。


s.t. : x+y\leq6
3x-2y\leq5
代码实现：
PSO参数设置
  1. # Work hard everyday!!!
  2. # Student:Sum3 Klaus
  3. # Time: 2022/3/13 16:31
  4. # Email: Sum3rebirth@gmail.com
  5. ”””
  6. pd.set_option(&#39;display.unicode.ambiguous_as_wide&#39;, True)  # 措置数据的列标题与数据无法对齐的情况
  7. pd.set_option(&#39;display.unicode.east_asian_width&#39;, True)  # 无法对齐主要是因为列标题是中文
  8. pd.set_option(&#39;display.max_columns&#39;, None)  # 显示所有列
  9. pd.set_option(&#39;display.max_rows&#39;, None)  # 显示所有行
  11. plt.rcParams[&#39;font.sans-serif&#39;] = [&#39;SimHei&#39;]  # 用来正常显示中文标签
  12. plt.rcParams[&#39;axes.unicode_minus&#39;] = False  # 用来正常显示负号
  14. import warnings
  15. warnings.filterwarnings(&#39;ignore&#39;)
  16. ”””
  18. import numpy as np
  19. from numpy import random
  20. import matplotlib.pyplot as plt
  21. import pandas as pd
  23. plt.rcParams[&#39;font.sans-serif&#39;] = [&#39;SimHei&#39;]
  24. plt.rcParams[&#39;axes.unicode_minus&#39;] = False
  25. # pd.set_option(&#39;display.max_rows&#39;, None)  # 显示所有行
  26. pd.set_option(&#39;display.max_columns&#39;, None)
  27. pd.set_option(&#39;display.unicode.ambiguous_as_wide&#39;, True)
  28. pd.set_option(&#39;display.unicode.east_asian_width&#39;, True)
  30. # PSO 参数
  31. w = 1
  32. w_min = 0.4
  33. w_max = 0.9
  34. c1 = 2
  35. c2 = 2
  36. r1 = None
  37. r2 = None
  38. dim = 2
  39. size = 100
  40. iter_max_num = 1000
  41. max_vel = 0.5
  42. fitness_value_lst = []
复制代码
定义适应度值计算函数
  1. def calc_fitness(X):
  2.     ”””计算粒子的适应度， shape=size * 2”””
  3.     A = 10
  4.     pi = np.pi
  5.     x = X[0]
  6.     y = X[1]
  8.     return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)
复制代码
定义约束惩罚项计算函数，以及权重计算函数
  1. def calc_e1(X):
  2.     ”””计算第一个约束的惩罚项”””
  3.     e_1 = X[0] + X[1] - 6
  4.     return max(0, e_1)
  7. def calc_e2(X):
  8.     ”””计算第一个约束的惩罚项”””
  9.     e_2 = 3 * X[0] - 2 * X[1] - 5
  10.     return max(0, e_2)
  13. def calc_Lj(e_1, e_2):
  14.     ”””按照每个粒子的约束惩罚项计算 Lj 权重值， e_1 e_2是列向量，暗示每个粒子的第一个与第二个约束的惩罚项值”””
  15.     if (e_1.sum() + e_2.sum()) <= 0:
  16.         # 防止分母为0
  17.         return 0, 0
  18.     else:
  19.         L_1 = e_1.sum() / (e_1.sum() + e_2.sum())
  20.         L_2 = e_1.sum() / (e_1.sum() + e_2.sum())
  21.         return L_1, L_2
复制代码
定义粒子群算法的速度更新函数
  1. def velocity_update(V, X, pbest, gbest, w):
  2.     ”””
  3.     按照速度迭代公式，更新粒子的速度
  4.     :param V: 粒子当前的速度矩阵，shape = 100*2
  5.     :param X: 粒子当前的位置矩阵，shape = 100*2
  6.     :param pbest:每个粒子历史最优位置，shape = 100*2
  7.     :param gbest:种群的历史最优位置， shape = 1*2
  8.     :return:
  9.     ”””
  10.     r1 = random.random((size, 1))
  11.     r2 = random.random((size, 1))
  12.     V = w * V + c1 * r1 * (pbest * X) + c2 * r2 * (gbest * X)
  13.     # 边界设置，防止越界
  14.     V[V < -max_vel] = -max_vel
  15.     V[V > max_vel] = max_vel
  16.     return V
复制代码
定义粒子群算法的位置更新函数
  1. def position_upodate(X, V):
  2.     ”””
  3.     更新粒子的位置
  4.     :param X: 粒子当前的位置矩阵，shape = 100*2
  5.     :param V: 粒子当前的速度矩阵，shape = 100*2
  6.     :return:
  7.     ”””
  8.     return X + V
复制代码
每个粒子的pbest更新函数
  1. # 10**（-7） = 0 计算机的数值精度位置
  2. zero = 10 ** (-7)
  5. def update_pbest(pbest, pbest_fitness, pbest_e, xi, xi_fitness, xi_e):
  6.     ”””
  7.     判断是否需要更新粒子的历史最优位置
  8.     :param pbest: 粒子i的历史最优位置
  9.     :param pbest_fitness: 粒子i的历史最优位置对应的方针函数值
  10.     :param pbest_e: 粒子i的历史最优位置对应的惩罚项约束
  11.     :param xi: 当前位置
  12.     :param xi_fitness:当前位置的方针函数值
  13.     :param xi_e: 当前位置的惩罚项值
  14.     :return:
  15.     ”””
  16.     # 择优法则1：如果pbest与xi都未违反约束，则拔取适应度值小的作为pbest
  17.     if pbest_e <= zero and xi_e <= zero:
  18.         if pbest_fitness <= xi_fitness:
  19.             return pbest, pbest_fitness, pbest_e
  20.         else:
  21.             return xi, xi_fitness, xi_e
  23.     # 择优法则2.1：如果当前位置违反约束，而历史最优位置没有违反，则选择历史最优
  24.     if pbest_e < zero and xi_e > zero:
  25.         return pbest, pbest_fitness, pbest_e
  27.     # 择优法则2.2：如果当前位置违反约束，而历史最优位置没有违反，则选择历史最优
  28.     if pbest_e > zero and xi_e < zero:
  29.         return xi, xi_fitness, xi_e
  31.     # 择优法则3：如果都违反约束，则拔取适应度值小的
  32.     if pbest_fitness < xi_fitness:
  33.         return pbest, pbest_fitness, pbest_e
  34.     else:
  35.         return xi, xi_fitness, xi_e
复制代码
种群gbest更新函数
  1. def update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e):
  2.     ”””
  3.     更新全局最优
  4.     :param gbest: 上一次迭代后的全局最优位置
  5.     :param gbest_fitness: 上一次迭代后的全局最优位置的适应度值
  6.     :param gbest_e: 上一次迭代后的全局最优位置的约束惩罚项
  7.     :param pbest: 当前迭代的种群的最优位置
  8.     :param pbest_fitness: 当前迭代的种群的最优位置的适应度值
  9.     :param pbest_e: 当前迭代的种群的最优位置的惩罚项值
  10.     :return:
  11.     ”””
  12.     # 先对种群寻求约束惩罚项未0的最优个体，如果每个个体的约束惩罚项都大于0，则寻找适应度值最小的个体
  13.     pbest_array = np.concatenate([pbest, pbest_fitness.reshape(-1, 1), pbest_e.reshape(-1, 1)], axis=1)
  14.     pbest_2 = pbest_array[pbest_array[:, -1] <= zero]  # 找出没有违反约束的粒子
  15.     if len(pbest_2) > 0:
  16.         pbest_2 = pbest_2[pbest_2[:, 2].argsort()]  # 按照适应度值排序
  17.     else:
  18.         pbest_2 = pbest_array[pbest_2[:, 2].argsort()]
  20.     # 当前迭代的最优个体
  21.     pbest_new, pbest_new_fitness, pbest_new_e = pbest_2[0, :2], pbest_2[0, 2], pbest_2[0, 3]
  23.     # 当前最优与全局最优斗劲
  24.     # 如果两者都没有约束
  25.     if gbest_e <= zero and pbest_new_e <= zero:
  26.         if gbest_fitness < pbest_new_fitness:
  27.             return gbest, gbest_fitness, gbest_e
  28.         else:
  29.             return pbest_new, pbest_new_fitness, pbest_new_e
  31.     # 此中只有有一者违反约束
  32.     if gbest_e <= zero and pbest_new_e > zero:
  33.         return gbest, gbest_fitness, gbest_e
  34.     if gbest_e > zero and pbest_new_e <= zero:
  35.         return pbest_new, pbest_new_fitness, pbest_new_e
  37.     # 如果二者都违反约束，则直接拔取适应度值最小的
  38.     if gbest_fitness < pbest_new_fitness:
  39.         return gbest, gbest_fitness, gbest_e
  40.     else:
  41.         return pbest_new, pbest_new_fitness, pbest_new_e
复制代码
主函数
  1. # 初始化一个矩阵info(shape=7) 记录迭代过程
  2. # 列属性
  3. # 1. 种群每个粒子的历史最优位置对应的适应度
  4. # 2. 历史最优位置对应的惩罚项
  5. # 3. 当前适应度值 （+ 惩罚项）
  6. # 4. 当前方针函数值
  7. # 5. 约束1的惩罚项
  8. # 6. 约束2的惩罚项
  9. # 7. 惩罚项之和
  10. info = np.zeros((size, 7))
  12. # 初始化种群各个粒子的位置
  13. X = random.uniform(-0.5, 0.5, size=(size, dim))
  15. # 初始化种群各个粒子的速度
  16. V = random.uniform(-0.5, 0.5, size=(size, dim))
  18. # 初始化种群各个粒子的历史最佳位置为当前位置
  19. pbest = X
  21. # 计算每个粒子的适应度
  22. for i in range(size):
  23.     info[i, 3] = calc_fitness(X[i])
  24.     info[i, 4] = calc_e1(X[i])
  25.     info[i, 5] = calc_e2(X[i])
  27.     # 计算惩罚项的权重，以及适应度值
  28.     L1, L2 = calc_Lj(info[i, 4], info[i, 5])
  29.     info[:, 6] = L1*info[:, 4] + L2*info[:, 5]
  30.     info[:, 2] = info[:, 3] + L1*info[:, 4] + L2*info[:, 5]
  32. # 历史最优
  33. info[:, 0] = info[:, 2]
  34. info[:, 1] = info[:, 6]
  36. # 全局最优
  37. gbest_i = info[:, 0].argmin()      # 全局最优的粒子的序号
  38. gbest = X[gbest_i]
  39. gbest_fitness = info[gbest_i, 0]
  40. gbest_e = info[gbest_i, 1]
  42. # 记录迭代过程的最优适应度值
  43. fitness_value_lst.append(gbest_fitness)
  45. # 迭代主函数
  48. def main(X, V, size, pbest, gbest, gbest_fitness, gbest_e, info, iter_max_num, w_min, w_max):
  49.     for j in range(iter_max_num):
  50.         w_iter = w_max - j*(w_max * w_min)/iter_max_num
  51.   
  52.         V = velocity_update(V, X, pbest=pbest, gbest=gbest, w=w_iter)
  54.         X = position_upodate(X, V)
  56.         for i in range(size):
  57.             info[i, 3] = calc_fitness(X[i])
  58.             info[i, 4] = calc_e1(X[i])
  59.             info[i, 5] = calc_e2(X[i])
  61.             L1, L2 = calc_Lj(info[i, 3], info[i, 3])
  62.             info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5]
  63.             info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5]
  65.         for i in range(size):
  66.             pbest_i = pbest[i]
  67.             pbest_i_fitness = info[i, 0]
  68.             pbest_i_e = info[i, 1]
  69.             x_i = X[i]
  70.             x_i_fitness = info[i, 2]
  71.             x_i_e = info[i, 6]
  73.             pbest_i, pbest_i_fitness, pbest_i_e = update_pbest(pbest_i, pbest_i_fitness, pbest_i_e,
  74.                                                                x_i, x_i_fitness, x_i_e)
  75.             pbest[i] = pbest_i
  76.             info[i, 0] = pbest_i_fitness
  77.             info[i, 1] = pbest_i_e
  79.         pbest_fitness = info[:, 2]
  80.         pbest_e = info[:, 6]
  81.         gbest, gbest_fitness, gbest_e = update_gbest(gbest, gbest_fitness, gbest_e,
  82.                                                      pbest, pbest_fitness, pbest_e)
  84.         fitness_value_lst.append(gbest_fitness)
  85.         columns = [&#39;0：历史最佳位置对应的适应度值&#39;, &#39;1：历史最优位置对应的惩罚项&#39;, &#39;2：当前适应度值&#39;,
  86.                    &#39;3：当前方针函数值&#39;, &#39;4：约束1的惩罚项&#39;, &#39;5：约束2的惩罚项&#39;, &#39;6:惩罚项的和&#39;]
  87.         info_pd = pd.DataFrame(info, columns=columns)
  88.         # print(info_pd)
复制代码


图-1.迭代过程



图-2.info_pd

欢迎各位老师同学批评斧正 ！
建议使用jupyter跑代码，pd的呈现会更美不雅观。
<hr/>

  • 《python最优化算法实战》

本帖子中包含更多资源

您需要 登录 才可以下载或查看，没有账号？立即注册

×
回复

使用道具 举报

返回列表 发新帖
懒得打字嘛，点击右侧快捷回复 【右侧内容，后台自定义】
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

小黑屋|手机版|Unity开发者联盟 ( 粤ICP备20003399号 )

GMT+8, 2024-11-17 02:35 , Processed in 0.143761 second(s), 28 queries .

Powered by Discuz! X3.5 Licensed

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表