Python pso粒子群算法的目标函数如何设置为迭代函数
这是迭代求得mfrc2，如何用pso算法求(mfrc2-3.5)^2最小时的H22输入量
PSO算法流程：
- 初始化粒子群；
随机设置各粒子的位置和速度，默认粒子的初始位置为粒子最优位置，并根据所有粒子最优位置，选取群体最优位置。 - 判断是否达到迭代次数；
若没有达到，则跳转到步骤3）。否则，直接输出结果。 - 更新所有粒子的位置和速度；
- 计算各粒子的适应度值。
# Import libs
import numpy as np
import random as rd
import matplotlib.pyplot as plt
# Constant definition
MIN_POS = [-5, -5] # Minimum position of the particle
MAX_POS = [5, 5] # Maximum position of the particle
MIN_SPD = [-0.5, -0.5] # Minimum speed of the particle
MAX_SPD = [1, 1] # Maximum speed of the particle
C1_MIN = 0
C1_MAX = 1.5
C2_MIN = 0
C2_MAX = 1.5
W_MAX = 1.4
W_MIN = 0
然后是PSO类
# Class definition
class PSO():
"""
PSO class
"""
def __init__(self,iters=100,pcount=50,pdim=2,mode='min'):
"""
PSO initialization
------------------
"""
self.w = None # Inertia factor
self.c1 = None # Learning factor
self.c2 = None # Learning factor
self.iters = iters # Number of iterations
self.pcount = pcount # Number of particles
self.pdim = pdim # Particle dimension
self.gbpos = np.array([0.0]*pdim) # Group optimal position
self.mode = mode # The mode of PSO
self.cur_pos = np.zeros((pcount, pdim)) # Current position of the particle
self.cur_spd = np.zeros((pcount, pdim)) # Current speed of the particle
self.bpos = np.zeros((pcount, pdim)) # The optimal position of the particle
self.trace = [] # Record the function value of the optimal solution
def init_particles(self):
"""
init_particles function
-----------------------
"""
# Generating particle swarm
for i in range(self.pcount):
for j in range(self.pdim):
self.cur_pos[i,j] = rd.uniform(MIN_POS[j], MAX_POS[j])
self.cur_spd[i,j] = rd.uniform(MIN_SPD[j], MAX_SPD[j])
self.bpos[i,j] = self.cur_pos[i,j]
# Initial group optimal position
for i in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.cur_pos[i]) < self.fitness(self.gbpos):
gbpos = self.cur_pos[i]
elif self.mode == 'max':
if self.fitness(self.cur_pos[i]) > self.fitness(self.gbpos):
gbpos = self.cur_pos[i]
def fitness(self, x):
"""
fitness function
----------------
Parameter:
x :
"""
# Objective function
fitval = 5*np.cos(x[0]*x[1])+x[0]*x[1]+x[1]**3 # min
# Retyrn value
return fitval
def adaptive(self, t, p, c1, c2, w):
"""
"""
#w = 0.95 #0.9-1.2
if t == 0:
c1 = 0
c2 = 0
w = 0.95
else:
if self.mode == 'min':
# c1
if self.fitness(self.cur_pos[p]) > self.fitness(self.bpos[p]):
c1 = C1_MIN + (t/self.iters)*C1_MAX + np.random.uniform(0,0.1)
elif self.fitness(self.cur_pos[p]) <= self.fitness(self.bpos[p]):
c1 = c1
# c2
if self.fitness(self.bpos[p]) > self.fitness(self.gbpos):
c2 = C2_MIN + (t/self.iters)*C2_MAX + np.random.uniform(0,0.1)
elif self.fitness(self.bpos[p]) <= self.fitness(self.gbpos):
c2 = c2
# w
#c1 = C1_MAX - (C1_MAX-C1_MIN)*(t/self.iters)
#c2 = C2_MIN + (C2_MAX-C2_MIN)*(t/self.iters)
w = W_MAX - (W_MAX-W_MIN)*(t/self.iters)
elif self.mode == 'max':
pass
return c1, c2, w
def update(self, t):
"""
update function
---------------
Note that :
1. Update particle position
2. Update particle speed
3. Update particle optimal position
4. Update group optimal position
"""
# Part1 : Traverse the particle swarm
for i in range(self.pcount):
# Dynamic parameters
self.c1, self.c2, self.w = self.adaptive(t,i,self.c1,self.c2,self.w)
# Calculate the speed after particle iteration
# Update particle speed
self.cur_spd[i] = self.w*self.cur_spd[i] \
+self.c1*rd.uniform(0,1)*(self.bpos[i]-self.cur_pos[i])\
+self.c2*rd.uniform(0,1)*(self.gbpos - self.cur_pos[i])
for n in range(self.pdim):
if self.cur_spd[i,n] > MAX_SPD[n]:
self.cur_spd[i,n] = MAX_SPD[n]
elif self.cur_spd[i,n] < MIN_SPD[n]:
self.cur_spd[i,n] = MIN_SPD[n]
# Calculate the position after particle iteration
# Update particle position
self.cur_pos[i] = self.cur_pos[i] + self.cur_spd[i]
for n in range(self.pdim):
if self.cur_pos[i,n] > MAX_POS[n]:
self.cur_pos[i,n] = MAX_POS[n]
elif self.cur_pos[i,n] < MIN_POS[n]:
self.cur_pos[i,n] = MIN_POS[n]
# Part2 : Update particle optimal position
for k in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.cur_pos[k]) < self.fitness(self.bpos[k]):
self.bpos[k] = self.cur_pos[k]
elif self.mode == 'max':
if self.fitness(self.cur_pos[k]) > self.fitness(self.bpos[k]):
self.bpos[k] = self.cur_pos[k]
# Part3 : Update group optimal position
for k in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.bpos[k]) < self.fitness(self.gbpos):
self.gbpos = self.bpos[k]
elif self.mode == 'max':
if self.fitness(self.bpos[k]) > self.fitness(self.gbpos):
self.gbpos = self.bpos[k]
def run(self):
"""
run function
-------------
"""
# Initialize the particle swarm
self.init_particles()
# Iteration
for t in range(self.iters):
# Update all particle information
self.update(t)
#
self.trace.append(self.fitness(self.gbpos))
然后是main
def main():
"""
main function
"""
for i in range(1):
pso = PSO(iters=100,pcount=50,pdim=2, mode='min')
pso.run()
#
print('='*40)
print('= Optimal solution:')
print('= x=', pso.gbpos[0])
print('= y=', pso.gbpos[1])
print('= Function value:')
print('= f(x,y)=', pso.fitness(pso.gbpos))
#print(pso.w)
print('='*40)
#
plt.plot(pso.trace, 'r')
title = 'MIN: ' + str(pso.fitness(pso.gbpos))
plt.title(title)
plt.xlabel("Number of iterations")
plt.ylabel("Function values")
plt.show()
#
input('= Press any key to exit...')
print('='*40)
exit()
if __name__ == "__main__":
main()
代码能复制贴出来吗，下面的迭代是干嘛，这样迭代，岂不是一个H22 会对应多个mfrc2 了？
这个应该需要把你的迭代公式当作一个适应度函数封装好，然后调用PSO模型。将测试函数那里换成你需要求(mfrc2-3.5)^2最小时的H22输入量时的函数
PSO模型如下：
'''
粒子群算法--------应用到单目标函数中
编写时间：2022.4.11
'''
import numpy as np
import matplotlib.pyplot as plt
import math
import random
from matplotlib.font_manager import FontProperties
font_set = FontProperties(fname=r"c:\windows\fonts\simsun.ttc", size=12)
# 测试函数
def test_function(x):
y = x * np.sin(x * math.pi * 10) + 2
return y
# 初始化参数
def init_param():
Np=100
D=1
c1 = 1.5
c2 = 2.5
w = 0.5
iterations=100
value_up_range=2
value_down_range=-1
return Np,D,c1,c2,w,iterations,value_up_range,value_down_range
# 定义存储矩阵
def EmptMatrix(Np, D,iterations):
individualBest_fitness=np.zeros(shape=(Np, 1))
individual_best=np.zeros(shape=(Np, D))
bestfitness = np.zeros(shape=(iterations, 1))
return individualBest_fitness,individual_best,bestfitness
# 初始化位置x和速度v
def init_x_v(Np,D,value_up_range,value_down_range,):
x=value_down_range+(value_up_range-value_down_range)*np.random.random([Np,D])
v=value_down_range+(value_up_range-value_down_range)*np.random.random([Np,D])
return x,v
# 计算各个粒子的适应值、个体最优值
def cal_fitness_individualBest(Np,x,individualBest_fitness,individual_best):
for i in range(Np):
individualBest_fitness[i,]=test_function(x[i,])
individual_best[i,]=x[i,]
return individualBest_fitness,individual_best
# 保存全局最优
def save_global_optimal(x,Np):
pg = x[0,]
for i in range(1, Np):
if test_function(x[i,]) < test_function(pg):
pg = x[i,]
return pg
# 进入主循环---迭代
def PSO(Np,D,c1,c2,w,iterations,bestfitness,value_up_range,value_down_range,x,v,individualBest_fitness,individual_best,pg):
for t in range(iterations):
# 次循环
for i in range(Np):
# 更新位置x,速度v
v[i,]=w*v[i,]+c1*np.random.random()*(individual_best[i,]-x[i,])+c2*np.random.random()*(pg-x[i,])
x[i,]=x[i,]+v[i,]
# 防止越界操作
for j in range(D):
if x[i,j]<value_down_range:
x[i,j] =value_down_range
if x[i,j]>value_up_range:
x[i,j]=value_up_range
# 个体最优选择比较
if test_function(x[i,])<individualBest_fitness[i,]:
individualBest_fitness[i,]=test_function(x[i,])
individual_best[i,]=x[i,]
# 全局最优选择比较
if individualBest_fitness[i,]<test_function(pg):
pg=individual_best[i,]
#保存每一代的最优值
bestfitness[t]=test_function(pg)
return pg,bestfitness
# 测试运行
Np,D,c1,c2,w,iterations,value_up_range,value_down_range=init_param()
individualBest_fitness,individual_best,bestfitness=EmptMatrix(Np, D,iterations)
x,v=init_x_v(Np,D,value_up_range,value_down_range,)
individualBest_fitness,individual_best=cal_fitness_individualBest(Np,x,individualBest_fitness,individual_best)
pg=save_global_optimal(x,Np)
pg,bestfitness=PSO(Np,D,c1,c2,w,iterations,bestfitness,value_up_range,value_down_range,x,v,individualBest_fitness,individual_best,pg)
# 输出结果、作图
print("函数取最小值时x等于：")
print(pg)
print('----------------')
print("最优值为：" )
print(bestfitness[iterations-1])
plt.plot(bestfitness,'r-.')
plt.ylabel("object fitness value")
plt.xlabel("迭代次数",fontproperties=font_set)
plt.title("the run result of function")
plt.show()
https://blog.csdn.net/brilliantZC/article/details/123846525
可以参考一下
第一步：对粒子群的随机位置和速度进行初始设定，同时设定迭代次数。
第二步：计算每个粒子的适应度值。
第三步：对每个粒子，将其适应度值与所经历的最好位置pbest i的适应度值进行比较，若较好，则将其作为当前的个体最优位置。
第四步：对每个粒子，将其适应度值与全局所经历的最好位置gbestg的适应度值进行比较，若较好，则将其作为当前的全局最优位置。
第五步：根据速度、位置公式对粒子的速度和位置进行优化，从而更新粒子位置。
第六步：如未达到结束条件（通常为最大循环数或最小误差要求），则返回第二步
import numpy as np
import random
class PSO_model:
def __init__(self,w,c1,c2,r1,r2,N,D,M):
self.w = w # 惯性权值
self.c1=c1
self.c2=c2
self.r1=r1
self.r2=r2
self.N=N # 初始化种群数量个数
self.D=D # 搜索空间维度
self.M=M # 迭代的最大次数
self.x=np.zeros((self.N,self.D)) #粒子的初始位置
self.v=np.zeros((self.N,self.D)) #粒子的初始速度
self.pbest=np.zeros((self.N,self.D)) #个体最优值初始化
self.gbest=np.zeros((1,self.D)) #种群最优值
self.p_fit=np.zeros(self.N)
self.fit=1e8 #初始化全局最优适应度
# 目标函数，也是适应度函数（求最小化问题）
def function(self,x):
A = 10
x1=x[0]
x2=x[1]
Z = 2 * A + x1 ** 2 - A * np.cos(2 * np.pi * x1) + x2 ** 2 - A * np.cos(2 * np.pi * x2)
return Z
# 初始化种群
def init_pop(self):
for i in range(self.N):
for j in range(self.D):
self.x[i][j] = random.random()
self.v[i][j] = random.random()
self.pbest[i] = self.x[i] # 初始化个体的最优值
aim=self.function(self.x[i]) # 计算个体的适应度值
self.p_fit[i]=aim # 初始化个体的最优位置
if aim < self.fit: # 对个体适应度进行比较，计算出最优的种群适应度
self.fit = aim
self.gbest = self.x[i]
# 更新粒子的位置与速度
def update(self):
for t in range(self.M): # 在迭代次数M内进行循环
for i in range(self.N): # 对所有种群进行一次循环
aim=self.function(self.x[i]) # 计算一次目标函数的适应度
if aim<self.p_fit[i]: # 比较适应度大小，将小的负值给个体最优
self.p_fit[i]=aim
self.pbest[i]=self.x[i]
if self.p_fit[i]<self.fit: # 如果是个体最优再将和全体最优进行对比
self.gbest=self.x[i]
self.fit = self.p_fit[i]
for i in range(self.N): # 更新粒子的速度和位置
self.v[i]=self.w*self.v[i]+self.c1*self.r1*(self.pbest[i]-self.x[i])+ self.c2*self.r2*(self.gbest-self.x[i])
self.x[i]=self.x[i]+self.v[i]
print("最优值：",self.fit,"位置为：",self.gbest)
if __name__ == '__main__':
# w,c1,c2,r1,r2,N,D,M参数初始化
w=random.random()
c1=c2=2#一般设置为2
r1=0.7
r2=0.5
N=30
D=2
M=200
pso_object=PSO_model(w,c1,c2,r1,r2,N,D,M)
pso_object.init_pop()
pso_object.update()
Class definition
class PSO():
"""
PSO class
"""
def __init__(self,iters=100,pcount=50,pdim=2,mode='min'):
"""
PSO initialization
------------------
"""
self.w = None # Inertia factor
self.c1 = None # Learning factor
self.c2 = None # Learning factor
self.iters = iters # Number of iterations
self.pcount = pcount # Number of particles
self.pdim = pdim # Particle dimension
self.gbpos = np.array([0.0]*pdim) # Group optimal position
self.mode = mode # The mode of PSO
self.cur_pos = np.zeros((pcount, pdim)) # Current position of the particle
self.cur_spd = np.zeros((pcount, pdim)) # Current speed of the particle
self.bpos = np.zeros((pcount, pdim)) # The optimal position of the particle
self.trace = [] # Record the function value of the optimal solution
def init_particles(self):
"""
init_particles function
-----------------------
"""
# Generating particle swarm
for i in range(self.pcount):
for j in range(self.pdim):
self.cur_pos[i,j] = rd.uniform(MIN_POS[j], MAX_POS[j])
self.cur_spd[i,j] = rd.uniform(MIN_SPD[j], MAX_SPD[j])
self.bpos[i,j] = self.cur_pos[i,j]
# Initial group optimal position
for i in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.cur_pos[i]) < self.fitness(self.gbpos):
gbpos = self.cur_pos[i]
elif self.mode == 'max':
if self.fitness(self.cur_pos[i]) > self.fitness(self.gbpos):
gbpos = self.cur_pos[i]
def fitness(self, x):
"""
fitness function
----------------
Parameter:
x :
"""
# Objective function
fitval = 5*np.cos(x[0]*x[1])+x[0]*x[1]+x[1]**3 # min
# Retyrn value
return fitval
def adaptive(self, t, p, c1, c2, w):
"""
"""
#w = 0.95 #0.9-1.2
if t == 0:
c1 = 0
c2 = 0
w = 0.95
else:
if self.mode == 'min':
# c1
if self.fitness(self.cur_pos[p]) > self.fitness(self.bpos[p]):
c1 = C1_MIN + (t/self.iters)*C1_MAX + np.random.uniform(0,0.1)
elif self.fitness(self.cur_pos[p]) <= self.fitness(self.bpos[p]):
c1 = c1
# c2
if self.fitness(self.bpos[p]) > self.fitness(self.gbpos):
c2 = C2_MIN + (t/self.iters)*C2_MAX + np.random.uniform(0,0.1)
elif self.fitness(self.bpos[p]) <= self.fitness(self.gbpos):
c2 = c2
# w
#c1 = C1_MAX - (C1_MAX-C1_MIN)*(t/self.iters)
#c2 = C2_MIN + (C2_MAX-C2_MIN)*(t/self.iters)
w = W_MAX - (W_MAX-W_MIN)*(t/self.iters)
elif self.mode == 'max':
pass
return c1, c2, w
def update(self, t):
"""
update function
---------------
Note that :
1. Update particle position
2. Update particle speed
3. Update particle optimal position
4. Update group optimal position
"""
# Part1 : Traverse the particle swarm
for i in range(self.pcount):
# Dynamic parameters
self.c1, self.c2, self.w = self.adaptive(t,i,self.c1,self.c2,self.w)
# Calculate the speed after particle iteration
# Update particle speed
self.cur_spd[i] = self.w*self.cur_spd[i] \
+self.c1*rd.uniform(0,1)*(self.bpos[i]-self.cur_pos[i])\
+self.c2*rd.uniform(0,1)*(self.gbpos - self.cur_pos[i])
for n in range(self.pdim):
if self.cur_spd[i,n] > MAX_SPD[n]:
self.cur_spd[i,n] = MAX_SPD[n]
elif self.cur_spd[i,n] < MIN_SPD[n]:
self.cur_spd[i,n] = MIN_SPD[n]
# Calculate the position after particle iteration
# Update particle position
self.cur_pos[i] = self.cur_pos[i] + self.cur_spd[i]
for n in range(self.pdim):
if self.cur_pos[i,n] > MAX_POS[n]:
self.cur_pos[i,n] = MAX_POS[n]
elif self.cur_pos[i,n] < MIN_POS[n]:
self.cur_pos[i,n] = MIN_POS[n]
# Part2 : Update particle optimal position
for k in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.cur_pos[k]) < self.fitness(self.bpos[k]):
self.bpos[k] = self.cur_pos[k]
elif self.mode == 'max':
if self.fitness(self.cur_pos[k]) > self.fitness(self.bpos[k]):
self.bpos[k] = self.cur_pos[k]
# Part3 : Update group optimal position
for k in range(self.pcount):
if self.mode == 'min':
if self.fitness(self.bpos[k]) < self.fitness(self.gbpos):
self.gbpos = self.bpos[k]
elif self.mode == 'max':
if self.fitness(self.bpos[k]) > self.fitness(self.gbpos):
self.gbpos = self.bpos[k]
def run(self):
"""
run function
-------------
"""
# Initialize the particle swarm
self.init_particles()
# Iteration
for t in range(self.iters):
# Update all particle information
self.update(t)
#
self.trace.append(self.fitness(self.gbpos))