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

huanhuanlala

想详细了解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最优化算法实战》