deffunction(x):#2004年考研题(我这里只求极小值):x^2-6xy+10^2-2yz-z^2+18=0,求z=z(x,y)的极值点和极值#此题答案极小值3,极小值点(9,3)#因为看答案知道z+y>0,所以只试z+y>0的情况
z=(18+x[0]**2-6*x[0]*x[1]+11*x[1]**2)**0.5-x[1]#化成x1,x2即y表示z的形式return z
实现
classParticle:def__init__(self, dim, max_opsition=100, max_velocity=0.05):
self.position =[random.uniform(-max_opsition, max_opsition)for i inrange(dim)]#粒子初始位置
self.velocity =[random.uniform(-max_velocity, max_velocity)for i inrange(dim)]# 粒子初始速度
self.particle_best =[0.0for i inrange(dim)]# 初始粒子最优解defset_position(self, i, value):
self.position[i]= value
defget_position(self):return self.position
defset_particle_best(self, value):
self.particle_best = value
defget_particle_best(self):return self.particle_best
defset_velocity(self, i, value):
self.velocity[i]= value
defget_velocity(self):return self.velocity
import random
particle_num =100#粒子数目
dim =2
particles =[Particle(dim=2)for i inrange(particle_num)]# 初始粒子
global_best =[0.0for i inrange(dim)]# 全局最优解defupdate(particles,global_best,iter_num,a,W,n1,n2,max_velocity,function):"""
状态更新
"""for j inrange(iter_num):for i inrange(len(particles)):# 计算每个粒子的最适应度,就是代个解看它的y是多少,如果y小于当前的最优解的值,就更新最优解为这个解if function(particles[i].get_position())< function(particles[i].get_particle_best()):
particles[i].set_particle_best(particles[i].get_position())# 更新全局最优解同理if function(particles[i].get_position())< function(global_best):
global_best = particles[i].get_position()for k inrange(dim):# 更新惯性权重
w = W[1]-(j+1)*(W[1]-W[0])/iter_num
# 套速度更新的公式
particles[i].set_velocity(k, w * particles[i].get_velocity()[k]+ \
n1 * random.random()*(a*particles[i].get_particle_best()[k]- particles[i].get_position()[k]) \
+ n2 * random.random()*(global_best[k]- particles[i].get_position()[k]))# 限制速度,别让它飞了if particles[i].get_velocity()[k]> max_velocity:
particles[i].get_velocity()[k]= max_velocity
elif particles[i].get_velocity()[k]<-max_velocity:
particles[i].get_velocity()[k]=-max_velocity
# 套位置更新的公式
particles[i].set_position(k,particles[i].get_position()[k]+ particles[i].get_velocity()[k])return global_best
a =1.#活化因子
W =[0,1]#惯性权重
n1 =1.#认知学习因子
n2 =2.#社会学习因子
max_velocity =0.05#速度限制
iter_num =2000#迭代次数
result=update(particles,global_best,iter_num,a,W,n1,n2,max_velocity,function)print(result)#极值点,答案约为(9,3)print(function(result))#极小值,答案约为3