赞
踩
最近在看清风的数学建模视频,其中有两道题:
# 第一步,绘制函数图像 import numpy as np import matplotlib.pyplot as plt def func(x): return 11*np.sin(x)+7*np.cos(5*x) x0 = np.linspace(-3,3,1000) y0 = func(x0) fig,ax = plt.subplots() ax.plot(x0,y0) ax.set_title('y = 11sin(x)+7cos(5x)') # 第二步,设置粒子群算法的参数 n = 10 # 粒子数量 narvs = 1 # 变量个数 c1 = 2 # 个体学习因子 c2 = 2 # 社会学习因子 w = 0.9 # 惯性权重 K = 50 # 迭代次数 vmax = (3-(-3))*0.2 # 粒子的最大速度 x_lb = -3 # x的下界 x_ub = 3 # x的上界 # 第三步,初始化粒子 x = x_lb + (x_ub-x_lb)*np.random.rand(n) v = -vmax + 2*vmax*np.random.rand(n) # 第四步,计算适应度 fit = func(x) # 计算每个粒子的适应度 pbest = x # 初始化这n个例子迄今为止找到的最佳位置 ind = np.argmax(fit) # 找到适应度最大的那个粒子的下标 gbest = x[ind] gbest_total = np.zeros(K) # 第五步,更新粒子速度和位置 for j in range(K): # 外层循环,共K次 v = w*v + c1*np.random.rand(n)*(pbest-x) + c2*np.random.rand(n)*(gbest-x) v[np.where(v<-vmax)] = -vmax # 速度小于-vmax的元素赋值为-vmax v[np.where(v>vmax)] = vmax # 速度大于vmax的元素赋值为vmax x = x + v # 更新第i个粒子的位置 x[np.where(x<x_lb)] = x_lb x[np.where(x>x_ub)] = x_ub # 第六步,重新计算适应度并找到最优粒子 fit = func(x) # 重新计算n个粒子的适应度 for k in range(n): # 更新第k个粒子迄今为止找到的最佳位置 if fit[k]>func(pbest[k]): pbest[k] = x[k] if np.max(fit)>func(gbest): # 更新所有粒子迄今找到的最佳位置 gbest = x[np.argmax(fit)] gbest_total[j] = func(gbest) ax.scatter(x,fit,c='r',marker='x') ax.scatter(gbest,func(gbest),c='r') ax.plot(gbest*np.ones(10),np.linspace(-20,func(gbest),10),'--') ax.plot(np.linspace(-3,gbest,10),func(gbest)*np.ones(10),'--') ax.set_ylim(-20,20) ax.autoscale(axis='x',tight=True) fig,ax = plt.subplots() ax.plot(np.arange(K),gbest_total) plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签 ax.set_title('算法找到的最大值随优化次数的变化曲线') ax.set_xlabel('循环次数') ax.autoscale(axis='x',tight=True) print('找到的最优解为:',func(gbest))
找到的最优解为: 17.492618154736594
红色叉号
表示最后一次循环结束后粒子所在的位置;红色圆圈(用橙色和绿色的线标识横纵坐标值)
表示算法寻到的最优解。# 第一步,绘制函数图像 import numpy as np import matplotlib.pyplot as plt import mpl_toolkits.mplot3d def func(x,y): return x**2+y**2-x*y-10*x-4*y+60 x0 = 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 = 0.9 # 惯性权重 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次 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]))
找到的最优解为: 8.031080661489156
图中红色叉号
表示最后一次循环结束后粒子所在的位置。
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。