参考文档链接:scikit-opt
本章继续Python的优化算法系列。
优化算法,尤其是启发式的仿生智能算法在最近很火,它适用于解决管理学,运筹学,统计学里面的一些优化问题。比如线性规划,整数规划,动态规划,非线性约束规划,甚至是超参数搜索等等方向的问题。
但是一般的优化算法还是matlab里面用的多,Python相关代码较少。博主在参考了很多文章的代码和模块之后,决定学习 scikit-opt 这个模块。这个优化算法模块对新手很友好,代码简洁,上手简单。而且代码和官方文档是中国人写的,还有很多案例,学起来就没什么压力...
缺点是包装的算法种类目前还不算多,只有七种:(差分进化算法、遗传算法、粒子群算法、模拟退火算法、蚁群算法、鱼群算法、免疫优化算法) /(其实已经够用了)
本次带来的是 数学建模里面经常使用的粒子群算法的使用演示。数学原理就不多说了
首先安装模块,在cmd里面或者anaconda prompt里面输入:
pip install scikit-opt
这个包很小,很快就能装好。
首先定义我们要求解的问题
- def demo_func(x):
- x1, x2, x3 = x
- return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2
调用PSO求解
- from sko.PSO import PSO
-
- pso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=150, lb=[0, -1, 0.5], ub=[1, 1, 1], w=0.8, c1=0.5, c2=0.5)
- pso.run()
- print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
画出y随着迭代次数的变化图
- import matplotlib.pyplot as plt
- plt.figure(figsize=(6,3))
- plt.plot(pso.gbest_y_hist)
- plt.show()
假如你的非线性约束是个圆内的面积 (x[0] - 1) ** 2 + (x[1] - 0) ** 2 - 0.5 ** 2<=0
这样写
- def demo_func(x):
- x1, x2,x3 = x
- return x1 ** 2 + (x2 - 0.05) ** 2 + x3 ** 2
- constraint_ueq = (lambda x: (x[0] - 1) ** 2 + (x[1] - 0) ** 2 - 0.5 ** 2 ,)
- pso = PSO(func=demo_func, n_dim=3, pop=40, max_iter=5000, lb=[-2, -2,-2], ub=[2, 2,2]
- , constraint_ueq=constraint_ueq)
- pso.run()
- print('best_x is ', pso.gbest_x, 'best_y is', pso.gbest_y)
(粒子群算法没有等式约束,都是不等式)
输入参数
输出参数
求解问题
- import numpy as np
- from sko.PSO import PSO
-
-
- def demo_func(x):
- x1, x2 = x
- return -20 * np.exp(-0.2 * np.sqrt(0.5 * (x1 ** 2 + x2 ** 2))) - np.exp(
- 0.5 * (np.cos(2 * np.pi * x1) + np.cos(2 * np.pi * x2))) + 20 + np.e
-
-
- constraint_ueq = (
- lambda x: (x[0] - 1) ** 2 + (x[1] - 0) ** 2 - 0.5 ** 2
- ,
- )
-
- max_iter = 50
- pso = PSO(func=demo_func, n_dim=2, pop=40, max_iter=max_iter, lb=[-2, -2], ub=[2, 2]
- , constraint_ueq=constraint_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 FuncAnimation
-
- record_value = pso.record_value
- X_list, V_list = record_value['X'], record_value['V']
-
- fig, ax = plt.subplots(1, 1)
- ax.set_title('title', loc='center')
- line = ax.plot([], [], 'b.')
-
- X_grid, Y_grid = np.meshgrid(np.linspace(-2.0, 2.0, 40), np.linspace(-2.0, 2.0, 40))
- Z_grid = demo_func((X_grid, Y_grid))
- ax.contour(X_grid, Y_grid, Z_grid, 30)
-
- ax.set_xlim(-2, 2)
- ax.set_ylim(-2, 2)
-
- t = np.linspace(0, 2 * np.pi, 40)
- ax.plot(0.5 * np.cos(t) + 1, 0.5 * np.sin(t), color='r')
-
- plt.ion()
- p = plt.show()
-
-
- def update_scatter(frame):
- i, j = frame // 10, frame % 10
- ax.set_title('iter = ' + str(i))
- X_tmp = X_list[i] + V_list[i] * j / 10.0
- plt.setp(line, 'xdata', X_tmp[:, 0], 'ydata', X_tmp[:, 1])
- return line
-
-
- ani = FuncAnimation(fig, update_scatter, blit=True, interval=25, frames=max_iter * 10)
- plt.show()
-
- ani.save('pso.gif', writer='pillow')
效果