智能优化算法学习-粒子群算法（Particle Swarm Optimization）-2解决约束优化问题(附python代码)
想详细了解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参数设置
# Work hard everyday!!!
# Student:Sum3 Klaus
# Time: 2022/3/13 16:31
# Email: Sum3rebirth@gmail.com
”””
pd.set_option('display.unicode.ambiguous_as_wide', True)# 措置数据的列标题与数据无法对齐的情况
pd.set_option('display.unicode.east_asian_width', True)# 无法对齐主要是因为列标题是中文
pd.set_option('display.max_columns', None)# 显示所有列
pd.set_option('display.max_rows', None)# 显示所有行
plt.rcParams['font.sans-serif'] = ['SimHei']# 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False# 用来正常显示负号
import warnings
warnings.filterwarnings('ignore')
”””
import numpy as np
from numpy import random
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
# pd.set_option('display.max_rows', None)# 显示所有行
pd.set_option('display.max_columns', None)
pd.set_option('display.unicode.ambiguous_as_wide', True)
pd.set_option('display.unicode.east_asian_width', True)
# PSO 参数
w = 1
w_min = 0.4
w_max = 0.9
c1 = 2
c2 = 2
r1 = None
r2 = None
dim = 2
size = 100
iter_max_num = 1000
max_vel = 0.5
fitness_value_lst = []定义适应度值计算函数
def calc_fitness(X):
”””计算粒子的适应度， shape=size * 2”””
A = 10
pi = np.pi
x = X
y = X
return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)定义约束惩罚项计算函数，以及权重计算函数
def calc_e1(X):
”””计算第一个约束的惩罚项”””
e_1 = X + X - 6
return max(0, e_1)
def calc_e2(X):
”””计算第一个约束的惩罚项”””
e_2 = 3 * X - 2 * X - 5
return max(0, e_2)
def calc_Lj(e_1, e_2):
”””按照每个粒子的约束惩罚项计算 Lj 权重值， e_1 e_2是列向量，暗示每个粒子的第一个与第二个约束的惩罚项值”””
if (e_1.sum() + e_2.sum()) <= 0:
# 防止分母为0
return 0, 0
else:
L_1 = e_1.sum() / (e_1.sum() + e_2.sum())
L_2 = e_1.sum() / (e_1.sum() + e_2.sum())
return L_1, L_2定义粒子群算法的速度更新函数
def velocity_update(V, X, pbest, gbest, w):
”””
按照速度迭代公式，更新粒子的速度
:param V: 粒子当前的速度矩阵，shape = 100*2
:param X: 粒子当前的位置矩阵，shape = 100*2
:param pbest:每个粒子历史最优位置，shape = 100*2
:param gbest:种群的历史最优位置， shape = 1*2
:return:
”””
r1 = random.random((size, 1))
r2 = random.random((size, 1))
V = w * V + c1 * r1 * (pbest * X) + c2 * r2 * (gbest * X)
# 边界设置，防止越界
V = -max_vel
V = max_vel
return V定义粒子群算法的位置更新函数
def position_upodate(X, V):
”””
更新粒子的位置
:param X: 粒子当前的位置矩阵，shape = 100*2
:param V: 粒子当前的速度矩阵，shape = 100*2
:return:
”””
return X + V每个粒子的pbest更新函数
# 10**（-7） = 0 计算机的数值精度位置
zero = 10 ** (-7)
def update_pbest(pbest, pbest_fitness, pbest_e, xi, xi_fitness, xi_e):
”””
判断是否需要更新粒子的历史最优位置
:param pbest: 粒子i的历史最优位置
:param pbest_fitness: 粒子i的历史最优位置对应的方针函数值
:param pbest_e: 粒子i的历史最优位置对应的惩罚项约束
:param xi: 当前位置
:param xi_fitness:当前位置的方针函数值
:param xi_e: 当前位置的惩罚项值
:return:
”””
# 择优法则1：如果pbest与xi都未违反约束，则拔取适应度值小的作为pbest
if pbest_e <= zero and xi_e <= zero:
if pbest_fitness <= xi_fitness:
return pbest, pbest_fitness, pbest_e
else:
return xi, xi_fitness, xi_e
# 择优法则2.1：如果当前位置违反约束，而历史最优位置没有违反，则选择历史最优
if pbest_e < zero and xi_e > zero:
return pbest, pbest_fitness, pbest_e
# 择优法则2.2：如果当前位置违反约束，而历史最优位置没有违反，则选择历史最优
if pbest_e > zero and xi_e < zero:
return xi, xi_fitness, xi_e
# 择优法则3：如果都违反约束，则拔取适应度值小的
if pbest_fitness < xi_fitness:
return pbest, pbest_fitness, pbest_e
else:
return xi, xi_fitness, xi_e种群gbest更新函数
def update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e):
”””
更新全局最优
:param gbest: 上一次迭代后的全局最优位置
:param gbest_fitness: 上一次迭代后的全局最优位置的适应度值
:param gbest_e: 上一次迭代后的全局最优位置的约束惩罚项
:param pbest: 当前迭代的种群的最优位置
:param pbest_fitness: 当前迭代的种群的最优位置的适应度值
:param pbest_e: 当前迭代的种群的最优位置的惩罚项值
:return:
”””
# 先对种群寻求约束惩罚项未0的最优个体，如果每个个体的约束惩罚项都大于0，则寻找适应度值最小的个体
pbest_array = np.concatenate(, axis=1)
pbest_2 = pbest_array <= zero]# 找出没有违反约束的粒子
if len(pbest_2) > 0:
pbest_2 = pbest_2.argsort()]# 按照适应度值排序
else:
pbest_2 = pbest_array.argsort()]
# 当前迭代的最优个体
pbest_new, pbest_new_fitness, pbest_new_e = pbest_2, pbest_2, pbest_2
# 当前最优与全局最优斗劲
# 如果两者都没有约束
if gbest_e <= zero and pbest_new_e <= zero:
if gbest_fitness < pbest_new_fitness:
return gbest, gbest_fitness, gbest_e
else:
return pbest_new, pbest_new_fitness, pbest_new_e
# 此中只有有一者违反约束
if gbest_e <= zero and pbest_new_e > zero:
return gbest, gbest_fitness, gbest_e
if gbest_e > zero and pbest_new_e <= zero:
return pbest_new, pbest_new_fitness, pbest_new_e
# 如果二者都违反约束，则直接拔取适应度值最小的
if gbest_fitness < pbest_new_fitness:
return gbest, gbest_fitness, gbest_e
else:
return pbest_new, pbest_new_fitness, pbest_new_e
主函数
# 初始化一个矩阵info(shape=7) 记录迭代过程
# 列属性
# 1. 种群每个粒子的历史最优位置对应的适应度
# 2. 历史最优位置对应的惩罚项
# 3. 当前适应度值 （+ 惩罚项）
# 4. 当前方针函数值
# 5. 约束1的惩罚项
# 6. 约束2的惩罚项
# 7. 惩罚项之和
info = np.zeros((size, 7))
# 初始化种群各个粒子的位置
X = random.uniform(-0.5, 0.5, size=(size, dim))
# 初始化种群各个粒子的速度
V = random.uniform(-0.5, 0.5, size=(size, dim))
# 初始化种群各个粒子的历史最佳位置为当前位置
pbest = X
# 计算每个粒子的适应度
for i in range(size):
info = calc_fitness(X)
info = calc_e1(X)
info = calc_e2(X)
# 计算惩罚项的权重，以及适应度值
L1, L2 = calc_Lj(info, info)
info[:, 6] = L1*info[:, 4] + L2*info[:, 5]
info[:, 2] = info[:, 3] + L1*info[:, 4] + L2*info[:, 5]
# 历史最优
info[:, 0] = info[:, 2]
info[:, 1] = info[:, 6]
# 全局最优
gbest_i = info[:, 0].argmin() # 全局最优的粒子的序号
gbest = X
gbest_fitness = info
gbest_e = info
# 记录迭代过程的最优适应度值
fitness_value_lst.append(gbest_fitness)
# 迭代主函数
def main(X, V, size, pbest, gbest, gbest_fitness, gbest_e, info, iter_max_num, w_min, w_max):
for j in range(iter_max_num):
w_iter = w_max - j*(w_max * w_min)/iter_max_num
V = velocity_update(V, X, pbest=pbest, gbest=gbest, w=w_iter)
X = position_upodate(X, V)
for i in range(size):
info = calc_fitness(X)
info = calc_e1(X)
info = calc_e2(X)
L1, L2 = calc_Lj(info, info)
info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5]
info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5]
for i in range(size):
pbest_i = pbest
pbest_i_fitness = info
pbest_i_e = info
x_i = X
x_i_fitness = info
x_i_e = info
pbest_i, pbest_i_fitness, pbest_i_e = update_pbest(pbest_i, pbest_i_fitness, pbest_i_e,
x_i, x_i_fitness, x_i_e)
pbest = pbest_i
info = pbest_i_fitness
info = pbest_i_e
pbest_fitness = info[:, 2]
pbest_e = info[:, 6]
gbest, gbest_fitness, gbest_e = update_gbest(gbest, gbest_fitness, gbest_e,
pbest, pbest_fitness, pbest_e)
fitness_value_lst.append(gbest_fitness)
columns = ['0：历史最佳位置对应的适应度值', '1：历史最优位置对应的惩罚项', '2：当前适应度值',
'3：当前方针函数值', '4：约束1的惩罚项', '5：约束2的惩罚项', '6:惩罚项的和']
info_pd = pd.DataFrame(info, columns=columns)
# print(info_pd)
图-1.迭代过程
图-2.info_pd
欢迎各位老师同学批评斧正 ！
建议使用jupyter跑代码，pd的呈现会更美不雅观。
<hr/>
[*]《python最优化算法实战》
页:
[1]