旧域名怎么做新网站,微信支付wordpress,做网站的域名,中国企业500强2018粒子群优化算法的实践
flyfish
粒子群优化算法#xff08;Particle Swarm Optimization#xff0c;PSO#xff09;或者粒子群算法 红叉的地方是理想之地#xff0c;这些粒子都想去#xff0c;总结8个字是信息共享#xff0c;个人决策。 上完图之后#xff0c;上代码Particle Swarm OptimizationPSO或者粒子群算法 红叉的地方是理想之地这些粒子都想去总结8个字是信息共享个人决策。 上完图之后上代码再解释
纯纯的PSO实现
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
np.set_printoptions(precision 4)
# 目标函数
def f(x,y):return (x-3.14)**2 (y-3.14)**2 np.sin(3*x3.14) np.cos(4*y-3.14)# 绘制三维函数
x, y np.array(np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100)))
z f(x, y)# 求出全局最小值
x_min np.around(x.ravel()[z.argmin()],4)
y_min np.around(y.ravel()[z.argmin()],4)
print(x_min:,x_min)
print(y_min:,y_min)
# 算法的超参数
c1 c2 0.1 #个体记忆 #集体记忆
w 0.8 #惯性权重 inertia weight constant.# 创建 particles
n_particles 20 # 20个粒子
np.random.seed(100)
X np.random.rand(2, n_particles) * 5
V np.random.randn(2, n_particles) * 0.1print(X:,X)
print(V:,V)# 初始化数据
pbest X # p personal cognitive
pbest_obj f(X[0], X[1])
gbest pbest[:, pbest_obj.argmin()] #g global social
gbest_obj pbest_obj.min()# 迭代一次粒子群优化
def update():global V, X, pbest, pbest_obj, gbest, gbest_obj# 更新参数r1, r2 np.random.rand(2)V w * V c1*r1*(pbest - X) c2*r2*(gbest.reshape(-1,1)-X)X X Vobj f(X[0], X[1])pbest[:, (pbest_obj obj)] X[:, (pbest_obj obj)]pbest_obj np.array([pbest_obj, obj]).min(axis0)gbest pbest[:, pbest_obj.argmin()]gbest_obj pbest_obj.min()# 等高线图
fig, ax plt.subplots(figsize(8,6))
fig.set_tight_layout(True)
img ax.imshow(z, extent[0, 5, 0, 5], originlower, cmapviridis, alpha0.5)
fig.colorbar(img, axax)
ax.plot([x_min], [y_min], markerx, markersize5, colorred)
contours ax.contour(x, y, z, 10, colorsblack, alpha0.4)
ax.clabel(contours, inlineTrue, fontsize8, fmt%.0f)
pbest_plot ax.scatter(pbest[0], pbest[1], markero, colorblack, alpha0.5)
p_plot ax.scatter(X[0], X[1], markero, colorblue, alpha0.5)
p_arrow ax.quiver(X[0], X[1], V[0], V[1], colorblue, width0.005, anglesxy, scale_unitsxy, scale1)
gbest_plot plt.scatter([gbest[0]], [gbest[1]], marker*, s100, colorblack, alpha0.4)
ax.set_xlim([0,5])
ax.set_ylim([0,5])# 粒子群算法的步骤算法更新和图形显示
def animate(i):title Iteration {:02d}.format(i)# 更新参数update()# 绘图ax.set_title(title)pbest_plot.set_offsets(pbest.T)p_plot.set_offsets(X.T)p_arrow.set_offsets(X.T)p_arrow.set_UVC(V[0], V[1])gbest_plot.set_offsets(gbest.reshape(1,-1))return ax, pbest_plot, p_plot, p_arrow, gbest_plotanim FuncAnimation(fig, animate, frameslist(range(1,20)), interval200, blitFalse, repeatTrue)
anim.save(PSO_1.gif, dpi120, writerimagemagick)print(PSO found best solution at f({}){}.format(gbest, ((gbest_obj))))
print(Global optimal at f({}){}.format([x_min,y_min], (f(x_min,y_min))))输出
x_min: 2.7273
y_min: 3.1313
X: [[2.717 1.3918 2.1226 4.2239 0.0236 0.6078 3.3537 4.1293 0.6835 2.87554.4566 1.046 0.9266 0.5419 1.0985 4.8931 4.0584 0.8597 4.0811 1.3704][2.1585 4.7001 4.0882 1.6806 0.8771 1.8642 0.0284 1.2621 3.9783 0.07632.9942 3.019 0.5257 1.9097 0.1824 4.4521 4.9046 0.2997 4.4527 2.8845]]
V: [[ 0.0731 0.1362 -0.0326 0.0056 0.0222 -0.1443 -0.0756 0.0816 0.075-0.0456 0.119 -0.1691 -0.1356 -0.1232 -0.0544 -0.0668 0.0007 -0.06130.13 -0.1733][-0.0983 0.0358 -0.1614 0.1471 -0.1188 -0.055 -0.094 -0.0828 0.01090.0508 -0.0862 0.1249 -0.008 -0.089 -0.0882 0.0019 0.0238 0.0014-0.1636 -0.1044]]
PSO found best solution at f([2.7157 3.1329])-1.7771745572809252
Global optimal at f([2.7273, 3.1313])-1.7760464968972247目标函数 f ( x , y ) ( x − 3.14 ) 2 ( y − 3.14 ) 2 sin ( 3 x 1.41 ) cos ( 4 y − 3.14 ) f(x,y) (x-3.14)^2 (y-3.14)^2 \sin(3x1.41) \cos(4y-3.14) f(x,y)(x−3.14)2(y−3.14)2sin(3x1.41)cos(4y−3.14)
那么我们如何才能找到这个函数的极小点呢当然我们可以采取穷举式搜索如果我们检查面上每个点的值我们就可以找到最小点。或者在面上随机找到一些样本点看看哪个样本点给出的值最低。。然而从形状上也可以注意到如果找到一个值较小的点则更容易在其附近找到更小值的点。
假设我们有粒子我们将粒子 i i i在迭代中的位置表示为 X i ( t ) X^i(t) Xi(t) 在上面的示例中我们将其作为坐标 X i ( t ) ( x i ( t ) , y i ( t ) ) . X^i(t) (x^i(t), y^i(t)). Xi(t)(xi(t),yi(t)).。 除了位置我们还有每个粒子的速度表示为 V i ( t ) ( v x i ( t ) , v y i ( t ) ) V^i(t)(v_x^i(t), v_y^i(t)) Vi(t)(vxi(t),vyi(t))。 在下一次迭代中每个粒子的位置将更新为 X i ( t 1 ) X i ( t ) V i ( t 1 ) X^i(t1) X^i(t)V^i(t1) Xi(t1)Xi(t)Vi(t1)
将大写X,拆成两个小写坐标下x,y表示如下 x i ( t 1 ) x i ( t ) v x i ( t 1 ) y i ( t 1 ) y i ( t ) v y i ( t 1 ) \begin{aligned} x^i(t1) x^i(t) v_x^i(t1) \\ y^i(t1) y^i(t) v_y^i(t1) \end{aligned} xi(t1)yi(t1)xi(t)vxi(t1)yi(t)vyi(t1)
同时速度也根据规则进行更新 V i ( t 1 ) w V i ( t ) c 1 r 1 ( p b e s t i – X i ( t ) ) c 2 r 2 ( g b e s t – X i ( t ) ) V^i(t1) w V^i(t) c_1r_1(pbest^i – X^i(t)) c_2r_2(gbest – X^i(t)) Vi(t1)wVi(t)c1r1(pbesti–Xi(t))c2r2(gbest–Xi(t)) 同样的公式还有如下 方式1 v i j ( t 1 ) w ∗ v i j ( t ) c p r 1 j ( t ) [ y i j ( t ) − x i j ( t ) ] c g r 2 j ( t ) [ y ^ j ( t ) − x i j ( t ) ] v_{ij}(t 1) w * v_{ij}(t) c_{p}r_{1j}(t)[y_{ij}(t) − x_{ij}(t)] c_{g}r_{2j}(t)[\hat{y}_{j}(t) − x_{ij}(t)] vij(t1)w∗vij(t)cpr1j(t)[yij(t)−xij(t)]cgr2j(t)[y^j(t)−xij(t)] 方式2 r1和r2是介于0和1之间的随机数,w,r1,r2是PSO算法的参数 X X X在上述中拆成了小写的 x , y x,y x,y,在代码中就是 X [ 0 ] , X [ 1 ] X[0], X[1] X[0],X[1] p b e s t i pbest^i pbesti的位置是给出粒子 i i i探索的最好 f ( X ) f(X) f(X)值Cognitivepersonal ,代码中的 p b e s t pbest pbest g b e s t gbest gbest是由群体中的所有粒子探索的Social(global)代码中的 g b e s t gbest gbest
请注意 p b e s t i pbest^i pbesti和 X i ( t ) ( x i ( t ) , y i ( t ) ) . X^i(t) (x^i(t), y^i(t)). Xi(t)(xi(t),yi(t)).是两个位置向量 X i ( t ) ( x i ( t ) , y i ( t ) ) X^i(t) (x^i(t), y^i(t)) Xi(t)(xi(t),yi(t)) p b e s t i ( x i ( t ) , y i ( t ) ) pbest^i(x^i(t), y^i(t)) pbesti(xi(t),yi(t)) 粒子 i i i最好的
差异 p b e s t i – X i ( t ) pbest^i – X^i(t) pbesti–Xi(t)是向量减法。 使用scikit-opt库实现
pip install scikit-opt该库封装了7种启发式算法差分进化算法、遗传算法、粒子群算法、模拟退火算法、蚁群算法、鱼群算法、免疫优化算法 源码地址
https://github.com/guofei9987/scikit-opt/输入
入参 默认值 意义
func - 目标函数
n_dim - 目标函数的维度
size_pop 50 种群规模
max_iter 200 最大迭代次数
lb None 每个参数的最小值
ub None 每个参数的最大值
w 0.8 惯性权重
c1 0.5 个体记忆
c2 0.5 集体记忆
constraint_ueq 空元组 不等式约束输出
pso.record_value 每一代的粒子位置、粒子速度、对应的函数值。pso.record_mode True 才开启记录
pso.gbest_y_hist 历史最优函数值
pso.best_y 最优函数值 迭代中使用的是 pso.gbest_x, pso.gbest_y
pso.best_x 最优函数值对应的输入值import numpy as np
from sko.PSO import PSOdef demo_func(X):x, y Xreturn (x-3.14)**2 (y-3.14)**2 np.sin(3*x3.14) np.cos(4*y-3.14)max_iter 30
# lb None 每个参数的最小值
# ub None 每个参数的最大值
pso PSO(funcdemo_func, n_dim2, pop40, max_itermax_iter, lb[-1, -1], ub[5, 5])pso.record_mode True
pso.run()
print(best_x is , pso.gbest_x, best_y is, pso.gbest_y)import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimationrecord_value pso.record_value
X_list, V_list record_value[X], record_value[V]fig, ax plt.subplots(1, 1)
ax.set_title(title, loccenter)
line ax.plot([], [], b.)X_grid, Y_grid np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100))Z_grid demo_func((X_grid, Y_grid))
ax.contour(X_grid, Y_grid, Z_grid, 10)ax.set_xlim([0,5])
ax.set_ylim([0,5])plt.ion()
p plt.show()def update_scatter(frame):i, j frame // 10, frame % 10ax.set_title(iter str(i))X_tmp X_list[i] V_list[i] * j / 10.0plt.setp(line, xdata, X_tmp[:, 0], ydata, X_tmp[:, 1])return lineani FuncAnimation(fig, update_scatter, blitTrue, interval25, framesmax_iter * 10)
plt.show()ani.save(PSO_2.gif, writerpillow)结果
best_x is [2.71374964 3.14188032] best_y is [-1.77777509]增加约束条件实现
在图上画个红色的紧箍,表示理想之地在红色圈内去那里吧
import numpy as np
from sko.PSO import PSOdef demo_func(X):x, y Xreturn (x-3.14)**2 (y-3.14)**2 np.sin(3*x3.14) np.cos(4*y-3.14)#非线性约束 (x[0] - 1) ** 2 (x[1] - 1) ** 2 - 10
constraint_ueq (lambda x: (x[0] - 1) ** 2 (x[1] - 1) ** 2-1,
)max_iter 30
# lb None 每个参数的最小值
# ub None 每个参数的最大值
pso PSO(funcdemo_func, n_dim2, pop40, max_itermax_iter, lb[-1, -1], ub[5, 5],constraint_ueqconstraint_ueq)pso.record_mode True
pso.run()
print(best_x is , pso.gbest_x, best_y is, pso.gbest_y)import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimationrecord_value pso.record_value
X_list, V_list record_value[X], record_value[V]fig, ax plt.subplots(1, 1)
ax.set_title(title, loccenter)
line ax.plot([], [], b.)X_grid, Y_grid np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100))Z_grid demo_func((X_grid, Y_grid))
ax.contour(X_grid, Y_grid, Z_grid, 10)ax.set_xlim([0,5])
ax.set_ylim([0,5])draw_circleplt.Circle((1.8, 1.5), 0.5,fillFalse,colorr)
ax.add_artist(draw_circle)plt.ion()
p plt.show()def update_scatter(frame):i, j frame // 10, frame % 10ax.set_title(iter str(i))X_tmp X_list[i] V_list[i] * j / 10.0plt.setp(line, xdata, X_tmp[:, 0], ydata, X_tmp[:, 1])return lineani FuncAnimation(fig, update_scatter, blitTrue, interval25, framesmax_iter * 10)
plt.show()ani.save(PSO_2.gif, writerpillow)
