一、粒子群算法
1.算法定义
粒子群算法,也称粒子群优化算法或鸟群觅食算法(Particle Swarm Optimization),缩写为 PSO, 是近年来由J. Kennedy和R. C. Eberhart等开发的一种新的进化算法(Evolutionary Algorithm - EA)。PSO 算法属于进化算法的一种,和模拟退火算法相似,它也是从随机解出发,通过迭代寻找最优解,它也是通过适应度来评价解的品质,但它比遗传算法规则更为简单,它没有遗传算法的“交叉”(Crossover) 和“变异”(Mutation) 操作,它通过追随当前搜索到的最优值来寻找全局最优。
2.算法优势
其迭代过程受自身经验、种群经验、速度惯性影响,这种算法以其实现容易、精度高、收敛快等优点引起了学术界的重视,并且在解决实际问题中展示了其优越性,在高维特征的收敛上有较好的表现。
代码实现(随机权重改进):
# 第一步,绘制函数图像
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3ddef func(x,y):
return x**2+y**2-x*y-10*x-4*y+60x0 = np.linspace(-15,15,100)
y0 = np.linspace(-15,15,100)
x0,y0 = np.meshgrid(x0,y0)
z0 = func(x0,y0)
fig = plt.figure(constrained_layout=True)
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x0,y0,z0,cmap=plt.cm.viridis,alpha=0.7)
#ax.plot_wireframe(x0,y0,z0) # 另一种绘图方式
ax.set_title('$y = x_1^2+x_2^2-x_1x_2-10x_1-4x_2+60$')# 第二步,设置粒子群算法的参数
n = 30 # 粒子数量
narvs = 2 # 变量个数
c1 = 2 # 个体学习因子
c2 = 2 # 社会学习因子
w_max = 0.9 # 惯性权重
w_min = 0.4
K = 40 # 迭代次数
vxmax = np.array([(15-(-15))*0.2,(15-(-15))*0.2]) # 粒子在x方向的最大速度
x_lb = np.array([-15,-15]) # x和y的下界
x_ub = np.array([15,15]) # x和y的上界# 第三步,初始化粒子
x = x_lb + (x_ub-x_lb)*np.random.rand(n,narvs)
v = -vxmax + 2*vxmax*np.random.rand(n,narvs)# 第四步,计算适应度
fit = func(x[:,0],x[:,1]) # 计算每个粒子的适应度
pbest = x # 初始化这n个例子迄今为止找到的最佳位置
ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标
gbest = x[ind,:]
gbest_total = np.zeros(K)# 第五步,更新粒子速度和位置
for j in range(K): # 外层循环,共K次
w = np.random.rand() # 随机数作为惯性权重
for p in range(n):
v[p,:] = w*v[p,:] + c1*np.random.rand(narvs)*(pbest[p,:]-x[p,:]) + c2*np.random.rand(narvs)*(gbest-x[p,:])
loc_v = np.where(v<-vxmax)
v[loc_v] = -vxmax[loc_v[1]] # 速度小于-vmax的元素赋值为-vmax
loc_v = np.where(v>vxmax)
v[loc_v] = vxmax[loc_v[1]] # 速度大于vmax的元素赋值为vmax
x = x + v # 更新第i个粒子的位置
loc_x = np.where(x<x_lb)
x[loc_x] = x_lb[loc_x[1]]
loc_x = np.where(x>x_ub)
x[loc_x] = x_ub[loc_x[1]]
# 第六步,重新计算适应度并找到最优粒子
fit = func(x[:,0],x[:,1]) # 重新计算n个粒子的适应度
for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置
if fit[k]<func(pbest[k,0],pbest[k,1]):
pbest[k,:] = x[k,:]
if np.min(fit)<func(gbest[0],gbest[1]): # 更新所有粒子迄今找到的最佳位置
gbest = x[np.argmin(fit),:]
gbest_total[j] = func(gbest[0],gbest[1])ax.scatter(x[:,0],x[:,1],fit,c='r',marker='x')
fig,ax = plt.subplots()
ax.plot(np.arange(K),gbest_total)
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False #用来显示负号
ax.set_title('算法找到的最小值随优化次数的变化曲线')
ax.set_xlabel('循环次数')
ax.autoscale(axis='x',tight=True)
print('找到的最优解为:',func(gbest[0],gbest[1]))
原文链接:改进的粒子群算法(PSO)算法Python代码——随机惯性权重_随即权重的pso-CSDN博客
二、遗传算法
1. 算法定义
遗传算法是一种模拟自然选择和遗传机制的优化算法,通过模拟基因的变异、交叉和选择等操作,逐代演化产生新的解,最终找到全局最优解。
2. 遗传算法的基本步骤
遗传算法的基本步骤包括:
- 初始化种群: 随机生成初始解的种群。
- 适应度评估: 计算每个个体的适应度,即解的优劣程度。
- 选择操作: 根据适应度选择个体,将适应度高的个体更有可能被选中。
- 交叉操作: 选中的个体进行基因交叉,产生新的个体。
- 变异操作: 对新个体进行基因变异,引入新的基因信息。
- 更新种群: 根据选择、交叉和变异等操作更新种群。
- 重复迭代: 重复进行选择、交叉、变异等操作,直到满足停止条件。
3. 个体的编码方法
在遗传算法中,个体的编码方式通常包括二进制编码、实数编码、排列编码等。选择适当的编码方式取决于具体问题的特点。
4.代码实现
import matplotlib.pyplot as plt
import numpy as np
import random
# (-1, 2)
# 初始化原始种群
def ori_popular(num):
popular = []
for i in range(num):
x = random.uniform(-1, 2) # 在此范围内生成一个随机浮点数
popular.append(x)
return popular
# 编码,也就是由表现型到基因型,性征到染色体
def encode(popular): # popular应该是float类型的列表
popular_gene = []
for i in range(0, len(popular)):
data = int((popular[i]-(-1)) / 3 * 2**18) # 染色体序列为18bit
bin_data = bin(data) # 整形转换成二进制是以字符串的形式存在的
for j in range(len(bin_data)-2, 18): # 序列长度不足补0
bin_data = bin_data[0:2] + '0' + bin_data[2:]
popular_gene.append(bin_data)
return popular_gene
# 解码,即适应度函数。通过基因,即染色体得到个体的适应度值
def decode(popular_gene):
fitness = []
for i in range(len(popular_gene)):
x = (int(popular_gene[i], 2) / 2**18) * 3 - 1
value = x * np.sin(10 * np.pi * x) + 2
fitness.append(value)
return fitness
# 选择and交叉。选择用轮牌赌,交叉概率为0.66
def choice_ex(popular_gene):
fitness = decode(popular_gene)
sum_fit_value = 0
for i in range(len(fitness)):
sum_fit_value += fitness[i]
# 各个个体被选择的概率
probability = []
for i in range(len(fitness)):
probability.append(fitness[i]/sum_fit_value)
# 概率分布
probability_sum = []
for i in range(len(fitness)):
if i == 0:
probability_sum.append(probability[i])
else:
probability_sum.append(probability_sum[i-1] + probability[i])
# 选择
popular_new = []
for i in range(int(len(fitness)/2)):
temp = []
for j in range(2):
rand = random.uniform(0, 1) # 在0-1之间随机一个浮点数
for k in range(len(fitness)):
if k == 0:
if rand < probability_sum[k]:
temp.append(popular_gene[k])
else:
if (rand > probability_sum[k-1]) and (rand < probability_sum[k]):
temp.append(popular_gene[k])
# 交叉,交叉率为0.66。
is_change = random.randint(0, 2)
if is_change:
temp_s = temp[0][9:15]
temp[0] = temp[0][0:9] + temp[1][9:15] + temp[0][15:]
temp[1] = temp[1][0:9] + temp_s + temp[1][15:]
popular_new.append(temp[0])
popular_new.append(temp[1])
return popular_new
# 变异.概率为0.05
def variation(popular_new):
for i in range(len(popular_new)):
is_variation = random.uniform(0, 1)
# print([len(k) for k in popular_new])
if is_variation < 0.02:
rand = random.randint(2, 19)
if popular_new[i][rand] == '0':
popular_new[i] = popular_new[i][0:rand] + '1' + popular_new[i][rand+1:]
else:
popular_new[i] = popular_new[i][0:rand] + '0' + popular_new[i][rand+1:]
return popular_new
if __name__ == '__main__': # alt+enter
# 初始化原始种群, 一百个个体
num = 100
ori_popular = ori_popular(num)
# 得到原始种群的基因
ori_popular_gene = encode(ori_popular) # 18位基因
new_popular_gene = ori_popular_gene
y = []
for i in range(1000): # 迭代次数。繁殖1000代
new_popular_gene = choice_ex(new_popular_gene) # 选择和交叉
new_popular_gene = variation(new_popular_gene) # 变异
# 取当代所有个体适应度平均值
new_fitness = decode(new_popular_gene)
sum_new_fitness = 0
for j in new_fitness:
sum_new_fitness += j
y.append(sum_new_fitness/len(new_fitness))
# 画图
x = np.linspace(0, 1000, 1000)
fig = plt.figure() # 相当于一个画板
axis = fig.add_subplot(111) # 坐标轴
axis.plot(x, y)
plt.show()
三、模拟退火算法
1.算法定义
用固体退火模拟组合优化问题,将内能E模拟为目标函数值f,温度T演化成控制参数t,即得到解组合优化问题的模拟退火算法:由初始解i和控制参数初值t开始,对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代,并逐步衰减t值,算法终止时的当前解即为所得近似最优解。
退火过程由冷却进度表(Cooling Schedule)控制,包括控制参数的初值t及其衰减因子Δt、每个t值时的迭代次数L和停止条件Tf。而温度的作用就是来计算转移概率P的。当温度每次下降后,转移概率也发生变化,因此在所有温度下迭代L次的结果也都是不相同的。在每个温度下迭代L次来寻找当前温度下的最优解,然后降低温度继续寻找,直到到达终止温度,即转移概率P接近于0。
2.退火过程中参数控制
(1)初始的温度T(0)应选的足够高,使的所有转移状态都被接受。初始温度越高,获得高质量的解的概率越大,耗费的时间越长。
(2)退火速率,即温度下降,最简单的下降方式是指数式下降:
T(n) = α \alphaαT(n) ,n =1,2,3,…
其中α \alphaα是小于1的正数,一般取值为0.8到0.99之间。使的对每一温度,有足够的转移尝试,指数式下降的收敛速度比较慢。
(3)终止温度
如果温度下降到终止温度或者达到用户设定的阈值,则退火完成。
import math
from random import random
import matplotlib.pyplot as pltdef func(x, y): #函数优化问题
res= 4*x**2-2.1*x**4+x**6/3+x*y-4*y**2+4*y**4
return res
#x为公式里的x1,y为公式里面的x2
class SA:
def __init__(self, func, iter=100, T0=100, Tf=0.01, alpha=0.99):
self.func = func
self.iter = iter #内循环迭代次数,即为L =100
self.alpha = alpha #降温系数,alpha=0.99
self.T0 = T0 #初始温度T0为100
self.Tf = Tf #温度终值Tf为0.01
self.T = T0 #当前温度
self.x = [random() * 11 -5 for i in range(iter)] #随机生成100个x的值
self.y = [random() * 11 -5 for i in range(iter)] #随机生成100个y的值
self.most_best =[]
"""
random()这个函数取0到1之间的小数
如果你要取0-10之间的整数(包括0和10)就写成 (int)random()*11就可以了,11乘以零点多的数最大是10点多,最小是0点多
该实例中x1和x2的绝对值不超过5(包含整数5和-5),(random() * 11 -5)的结果是-6到6之间的任意值(不包括-6和6)
(random() * 10 -5)的结果是-5到5之间的任意值(不包括-5和5),所有先乘以11,取-6到6之间的值,产生新解过程中,用一个if条件语句把-5到5之间(包括整数5和-5)的筛选出来。
"""
self.history = {'f': [], 'T': []}def generate_new(self, x, y): #扰动产生新解的过程
while True:
x_new = x + self.T * (random() - random())
y_new = y + self.T * (random() - random())
if (-5 <= x_new <= 5) & (-5 <= y_new <= 5):
break #重复得到新解,直到产生的新解满足约束条件
return x_new, y_newdef Metrospolis(self, f, f_new): #Metropolis准则
if f_new <= f:
return 1
else:
p = math.exp((f - f_new) / self.T)
if random() < p:
return 1
else:
return 0def best(self): #获取最优目标函数值
f_list = [] #f_list数组保存每次迭代之后的值
for i in range(self.iter):
f = self.func(self.x[i], self.y[i])
f_list.append(f)
f_best = min(f_list)
idx = f_list.index(f_best)
return f_best, idx #f_best,idx分别为在该温度下,迭代L次之后目标函数的最优解和最优解的下标def run(self):
count = 0
#外循环迭代,当前温度小于终止温度的阈值
while self.T > self.Tf:
#内循环迭代100次
for i in range(self.iter):
f = self.func(self.x[i], self.y[i]) #f为迭代一次后的值
x_new, y_new = self.generate_new(self.x[i], self.y[i]) #产生新解
f_new = self.func(x_new, y_new) #产生新值
if self.Metrospolis(f, f_new): #判断是否接受新值
self.x[i] = x_new #如果接受新值,则把新值的x,y存入x数组和y数组
self.y[i] = y_new
# 迭代L次记录在该温度下最优解
ft, _ = self.best()
self.history['f'].append(ft)
self.history['T'].append(self.T)
#温度按照一定的比例下降(冷却)
self.T = self.T * self.alpha
count += 1
# 得到最优解
f_best, idx = self.best()
print(f"F={f_best}, x={self.x[idx]}, y={self.y[idx]}")sa = SA(func)
sa.run()plt.plot(sa.history['T'], sa.history['f'])
plt.title('SA')
plt.xlabel('T')
plt.ylabel('f')
plt.gca().invert_xaxis()
plt.show()