粒子群算法解决非线性问题
引入
上次我们介绍了粒子群算法的各种改进，以及matlab软件自带的更强大的粒子群算法，解决的问题都是连续的，无约束的；那么我们能解决有约束的，非线性问题吗？
当然可以，不过在此之前，我们需要搞清实现的思路。
解决非线性问题的两种思路
- 直接在更新新的个体位置之前加入约束条件，满足约束条件更新个体位置，不满足约束条件，不更新个体位置
- 加入惩罚项，引入罚函数。在适应度函数上构造惩罚项，对违反约束的情况进行惩罚，将有约束的优化问题转化为无约束的优化问题。
罚函数解决非线性约束问题理论
在解决非线性约束问题中最常见的手段是引入罚函数。
什么是罚函数？
罚函数是指在求解最优化问题(无线性约束优化及非线性约束优化)时，在原有目标函数中加上一个障碍函数，而得到一个增广目标函数，罚函数的功能是对非可行点或企图穿越边界而逃离可行域的点赋予一个极大的值，即将有约束最优化问题转化为求解无约束最优化问题。
简单来说，就是对对每一个粒子进行判断，判断该粒子的取值是否满足约束条件，如果满足，则进行正常的适应度计算，如果不满足，则将该粒子的适应度直接赋值为一个很大的数。这个很大的数我们叫做惩罚。
约束条件一般由等式约束和不等式约束组成。我们一般将等式约束转化为不等式约束处理。
加入惩罚项之后，我们新的目标函数变为：
其中\(σ\)是惩罚因子，一般取$σ=t√t $，P(x)是整体惩罚项，P(x)的计算方法如下：
1.对于不等式\(g_i(x)<=0\),惩罚项
即如果\(g_i(x)<=0\),则\(e_i(x)=0\),否则\(e_i(x)=-g_i(x)\)
2.对于等式约束,一个简单的方法设定一个等式约束容忍度值\(ε\),新的不等式约束为:
因此等式约束的惩罚项为:
整体惩罚性P(x)是各个约束惩罚项的和，即：
由于约束条件之间的差异，某些约束条件可能对个体违反程度起着决定性的作用，为了平衡这种差异，对惩罚项做归一化处理:
最后得到的粒子xi的整体惩罚项P(x)的计算公式：
在粒子群算法中，每一步迭代都会更新Pbest和Gbest，虽然可以将有约束问题转换为无约束问题进行迭代求解，但是问题的解xi依然存在不满足约束条件的情况，因此需要编制一些规则来比较两个粒子的优劣，规则如下：
- 如果两个粒子\(xi\)和\(xj\)都可行，则比较其适应度函数\(f(x_i)\)和\(f(x_j)\)，值小的粒子为优。
- 当两个粒子\(xi\)和\(xj\)都不可行，则比较惩罚项\(P(x_i)\)和\(P(x_j）\)，违背约束程度小的粒子更优。
- 当粒子\(x_i\)可行而粒子\(x_j\)不可行，选可行解。
对于粒子的上下限约束 可以体现在位置更新函数里，不必加惩罚项。 具体思路就是遍历每一个粒子的位置，如果超除上下限，位置则更改为上下限中的任何一个位置.
罚函数解决非线性约束问题实现
例子
初始化参数
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
# PSO的参数
w = 1 # 惯性因子，一般取1
c1 = 2 # 学习因子，一般取2
c2 = 2 #
r1 = None # 为两个（0,1）之间的随机数
r2 = None
dim = 2 # 维度的维度#对应2个参数x,y
size = 100 # 种群大小，即种群中小鸟的个数
iter_num = 1000 # 算法最大迭代次数
max_vel = 0.5 # 限制粒子的最大速度为0.5
fitneess_value_list = [] # 记录每次迭代过程中的种群适应度值变化
计算适应度函数和计算约束惩罚项函数
def calc_f(X):
"""计算粒子的的适应度值，也就是目标函数值，X 的维度是 size * 2 """
A = 10
pi = np.pi
x = X[0]
y = X[1]
return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)
def calc_e1(X):
"""计算第一个约束的惩罚项"""
e = X[0] + X[1] - 6
return max(0, e)
def calc_e2(X):
"""计算第二个约束的惩罚项"""
e = 3 * X[0] - 2 * X[1] - 5
return max(0, e)
def calc_Lj(e1, e2):
"""根据每个粒子的约束惩罚项计算Lj权重值，e1, e2列向量，表示每个粒子的第1个第2个约束的惩罚项值"""
# 注意防止分母为零的情况
if (e1.sum() + e2.sum()) <= 0:
return 0, 0
else:
L1 = e1.sum() / (e1.sum() + e2.sum())
L2 = e2.sum() / (e1.sum() + e2.sum())
return L1, L2
定义粒子群算法的速度更新函数，位置更新函数
def velocity_update(V, X, pbest, gbest):
"""
根据速度更新公式更新每个粒子的速度
种群size=20
:param V: 粒子当前的速度矩阵，20*2 的矩阵
:param X: 粒子当前的位置矩阵，20*2 的矩阵
:param pbest: 每个粒子历史最优位置，20*2 的矩阵
:param gbest: 种群历史最优位置，1*2 的矩阵
"""
r1 = np.random.random((size, 1))
r2 = np.random.random((size, 1))
V = w * V + c1 * r1 * (pbest - X) + c2 * r2 * (gbest - X) # 直接对照公式写就好了
# 防止越界处理
V[V < -max_vel] = -max_vel
V[V > max_vel] = max_vel
return V
def position_update(X, V):
"""
根据公式更新粒子的位置
:param X: 粒子当前的位置矩阵，维度是 20*2
:param V: 粒子当前的速度举着，维度是 20*2
"""
X=X+V#更新位置
size=np.shape(X)[0]#种群大小
for i in range(size):#遍历每一个例子
if X[i][0]<=1 or X[i][0]>=2:#x的上下限约束
X[i][0]=np.random.uniform(1,2,1)[0]#则在1到2随机生成一个数
if X[i][1] <= -1 or X[i][0] >= 0:#y的上下限约束
X[i][1] = np.random.uniform(-1, 0, 1)[0] # 则在-1到0随机生成一个数
return X
每个粒子历史最优位置更优函数，以及整个群体历史最优位置更新函数
def update_pbest(pbest, pbest_fitness, pbest_e, xi, xi_fitness, xi_e):
"""
判断是否需要更新粒子的历史最优位置
:param pbest: 历史最优位置
:param pbest_fitness: 历史最优位置对应的适应度值
:param pbest_e: 历史最优位置对应的约束惩罚项
:param xi: 当前位置
:param xi_fitness: 当前位置的适应度函数值
:param xi_e: 当前位置的约束惩罚项
:return:
"""
# 下面的 0.0000001 是考虑到计算机的数值精度位置，值等同于0
# 规则1，如果 pbest 和 xi 都没有违反约束，则取适应度小的
if pbest_e <= 0.0000001 and xi_e <= 0.0000001:
if pbest_fitness <= xi_fitness:
return pbest, pbest_fitness, pbest_e
else:
return xi, xi_fitness, xi_e
# 规则2，如果当前位置违反约束而历史最优没有违反约束，则取历史最优
if pbest_e < 0.0000001 and xi_e >= 0.0000001:
return pbest, pbest_fitness, pbest_e
# 规则3，如果历史位置违反约束而当前位置没有违反约束，则取当前位置
if pbest_e >= 0.0000001 and xi_e < 0.0000001:
return xi, xi_fitness, xi_e
# 规则4，如果两个都违反约束，则取适应度值小的
if pbest_fitness <= xi_fitness:
return pbest, pbest_fitness, pbest_e
else:
return xi, xi_fitness, xi_e
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，就找适应度最小的个体
pbest2 = np.concatenate([pbest, pbest_fitness.reshape(-1, 1), pbest_e.reshape(-1, 1)], axis=1) # 将几个矩阵拼接成矩阵 ，4维矩阵（x,y,fitness,e）
pbest2_1 = pbest2[pbest2[:, -1] <= 0.0000001] # 找出没有违反约束的个体
if len(pbest2_1) > 0:
pbest2_1 = pbest2_1[pbest2_1[:, 2].argsort()] # 根据适应度值排序
else:
pbest2_1 = pbest2[pbest2[:, 2].argsort()] # 如果所有个体都违反约束，直接找出适应度值最小的
# 当前迭代的最优个体
pbesti, pbesti_fitness, pbesti_e = pbest2_1[0, :2], pbest2_1[0, 2], pbest2_1[0, 3]
# 当前最优和全局最优比较
# 如果两者都没有约束
if gbest_e <= 0.0000001 and pbesti_e <= 0.0000001:
if gbest_fitness < pbesti_fitness:
return gbest, gbest_fitness, gbest_e
else:
return pbesti, pbesti_fitness, pbesti_e
# 有一个违反约束而另一个没有违反约束
if gbest_e <= 0.0000001 and pbesti_e > 0.0000001:
return gbest, gbest_fitness, gbest_e
if gbest_e > 0.0000001 and pbesti_e <= 0.0000001:
return pbesti, pbesti_fitness, pbesti_e
# 如果都违反约束，直接取适应度小的
if gbest_fitness < pbesti_fitness:
return gbest, gbest_fitness, gbest_e
else:
return pbesti, pbesti_fitness, pbesti_e
PSO
# 初始化一个矩阵 info, 记录：
# 0、种群每个粒子的历史最优位置对应的适应度，
# 1、历史最优位置对应的惩罚项，
# 2、当前适应度，
# 3、当前目标函数值，
# 4、约束1惩罚项，
# 5、约束2惩罚项，
# 6、惩罚项的和
# 所以列的维度是7
info = np.zeros((size, 7))
# 初始化种群的各个粒子的位置
# 用一个 20*2 的矩阵表示种群，每行表示一个粒子
X = np.random.uniform(-5, 5, size=(size, dim))
# 初始化种群的各个粒子的速度
V = np.random.uniform(-0.5, 0.5, size=(size, dim))
# 初始化粒子历史最优位置为当当前位置
pbest = X
# 计算每个粒子的适应度
for i in range(size):
info[i, 3] = calc_f(X[i]) # 目标函数值
info[i, 4] = calc_e1(X[i]) # 第一个约束的惩罚项
info[i, 5] = calc_e2(X[i]) # 第二个约束的惩罚项
# 计算惩罚项的权重，及适应度值
L1, L2 = calc_Lj(info[i, 4], info[i, 5])
info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5] # 适应度值
info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5] # 惩罚项的加权求和
# 历史最优
info[:, 0] = info[:, 2] # 粒子的历史最优位置对应的适应度值
info[:, 1] = info[:, 6] # 粒子的历史最优位置对应的惩罚项值
# 全局最优
gbest_i = info[:, 0].argmin() # 全局最优对应的粒子编号
gbest = X[gbest_i] # 全局最优粒子的位置
gbest_fitness = info[gbest_i, 0] # 全局最优位置对应的适应度值
gbest_e = info[gbest_i, 1] # 全局最优位置对应的惩罚项
# 记录迭代过程的最优适应度值
fitneess_value_list.append(gbest_fitness)
# 接下来开始迭代
for j in range(iter_num):
# 更新速度
V = velocity_update(V, X, pbest=pbest, gbest=gbest)
# 更新位置
X = position_update(X, V)
# 计算每个粒子的目标函数和约束惩罚项
for i in range(size):
info[i, 3] = calc_f(X[i]) # 目标函数值
info[i, 4] = calc_e1(X[i]) # 第一个约束的惩罚项
info[i, 5] = calc_e2(X[i]) # 第二个约束的惩罚项
# 计算惩罚项的权重，及适应度值
L1, L2 = calc_Lj(info[i, 4], info[i, 5])
info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5] # 适应度值
info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5] # 惩罚项的加权求和
# 更新历史最优位置
for i in range(size):
pbesti = pbest[i]
pbest_fitness = info[i, 0]
pbest_e = info[i, 1]
xi = X[i]
xi_fitness = info[i, 2]
xi_e = info[i, 6]
# 计算更新个体历史最优
pbesti, pbest_fitness, pbest_e = \
update_pbest(pbesti, pbest_fitness, pbest_e, xi, xi_fitness, xi_e)
pbest[i] = pbesti
info[i, 0] = pbest_fitness
info[i, 1] = pbest_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)
# 记录当前迭代全局之硬度
fitneess_value_list.append(gbest_fitness)
# 最后绘制适应度值曲线
print('迭代最优结果是：%.5f' % calc_f(gbest))
print('迭代最优变量是：x=%.5f, y=%.5f' % (gbest[0], gbest[1]))
print('迭代约束惩罚项是：', gbest_e)
#迭代最优结果是：1.00491
#迭代最优变量是：x=1.00167, y=-0.00226
#迭代约束惩罚项是： 0.0
# 从结果看，有多个不同的解的目标函数值是相同的，多测试几次就发现了
# 绘图
plt.plot(fitneess_value_list[: 30], color='r')
plt.title('迭代过程')
plt.show()
本文来自：粒子群算法求解带约束优化问题 源码实现_总裁余（余登武）博客-CSDN博客_带约束粒子群
可参考的其他相关博客：