python：scipy.optimize.minimize(method=’Nelder-Mead’)

简介

优化问题是工程实践中经常遇到的一种问题。简单讲，就是搜索优化出一组自变量参数，使得目标函数达到极小值(极大值)。

如何搜索出这组参数呢：这就是优化算法做的事情。不同的优化算法适用于不同的优化问题。

本文简要介绍在python种NM算法来解决局部优化问题。

接口

使用总结：

• NM是局域优化，即最终的优化结果会掉到一个局域最小值附近，故对初值敏感
• 有两种方式设置初值，1.设置x0；2.设置options中的initial_simplex。手动设置initial_simplex更好，可灵活控制收敛速度。
• 终止条件有：maxiter, maxfev, xatol, fatol。判断条件是：先判断是否达到maxiter, maxfev，若达到则终止；若未达到，则再判断 xatol, fatol是否同时满足条件，若满足，则终止。其源代码逻辑如下：

  while (fcalls[0] < maxfun and iterations < maxiter):
       if (numpy.max(numpy.ravel(numpy.abs(sim[1:] - sim[0]))) <= xatol and
               numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= fatol):
           break

• maxiter, maxfev默认值为N*200。比如自变量个数为2，则默认值为400；xatol, fatol默认值为1e-4，一般根据实际问题手动设置。越小，则收敛越快。
• 若不同参数的取值范围不同，则应对其进行归一化处理

实例1：Himmelblau函数

关于Himmelblau函数

$f(x,y)=(x^2+y-11{)}^2+(x+y^2-7{)}^2$
有四个局域最小值: [3,2], [-2.805,3.131], [-3.779,-3.283], [3.584,-1.848]。

def f(x,y):
    return (x**2+y-11)**2+(x+y**2-7)**2

x = np.linspace(-7,7,50)
y = np.linspace(-7,7,50)

z = np.zeros((50,50))
for i,a in enumerate(x):
    for j,b in enumerate(y):      
        z[i,j] = f(a,b)
xx,yy = np.meshgrid(x,y)
fig, ax = plt.subplots()
c = ax.pcolormesh(xx,yy,z.T,cmap='jet')
fig.colorbar(c, ax=ax)

优化流程

实例代码

• 注意，为了绘制出迭代过程，在cost_function里加了一个xs.append(x)。实际问题中可根据需求看是否需要。本实例中将迭代过程中的参数绘制到二维扫描图上，可形象的演示NM算法的搜索过程。

xs = []
def cost_function(x):
    xs.append(x)
    return (x[0]**2+x[1]-11)**2+(x[0]+x[1]**2-7)**2

x_center = np.array([0,0])
step = 0.5
x0 = np.vstack((x_center, x_center+np.diag((step,step))))
xatol,fatol = 1e-3,1e-3
res = minimize(cost_function,x_center,method='nelder-mead',options={'initial_simplex':x0,'xatol': xatol,'fatol': fatol, 'disp': True})
print(res.x)
x_ = np.array([x[0] for x in xs])
y_ = np.array([x[1] for x in xs])
ys = np.array([f(x[0],x[1]) for x in xs])

fig, ax = plt.subplots()
c = ax.pcolormesh(xx,yy,z.T,cmap='jet')
fig.colorbar(c, ax=ax)
ax.plot(x_,y_,'r',x_[0],y_[0],'go',x_[-1],y_[-1],'y+',markersize=6)

fig2 = plt.figure()
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax3 = plt.subplot(212)
ax1.plot(x_)
ax1.set_title('x')
ax2.plot(y_)
ax2.set_title('y')
ax3.plot(ys)

结果及分析：

• 迭代了35次，最终达到xatol,fatol的终止条件。
• 初值流形是在[0, 0]，NM是局域优化，故最终搜索到的最优值大约是[3.00007588 1.99958898]，符合理论中的一个局域最小值[3, 2]。

Optimization terminated successfully.
         Current function value: 0.000002
         Iterations: 35
         Function evaluations: 68
[3.00021471 1.99974856]

对初值敏感的局域优化

NM算法时局域优化算法，所以对初值敏感。下面分别选择四个初值：[0,0], [-1, 0],[0,-1],[-1,1] 进行优化。

代码同上，只修改初值，比如：x_center = np.array([-1,0])

结果：符合预期。
初值为[0,0]，搜索到[3,2]；
初值为[-1,0]，搜索到[-2.805,3.131]；
初值为[0,-1]，搜索到[3.584,-1.848]；
初值为[-1,-1]，搜索到[-3.779,-3.283]；

归一化问题

在实际问题中，有时不同的自变量参数的量纲并不同，即取值范围和搜索步长不同。此时应该进行归一化处理。

稍微修改函数，模拟成一个需要归一化的优化问题。

实例代码

def f2(x,y):
    return (x**2+y/1000-11)**2+(x+(y/1000)**2-7)**2

x = np.linspace(-7,7,50)
y = np.linspace(-7000,7000,50)

z = np.zeros((50,50))
for i,a in enumerate(x):
    for j,b in enumerate(y):      
        z[i,j] = f2(a,b)
xx,yy = np.meshgrid(x,y)

xs = []
def cost_function2(x):
    xs.append(x)
    x_norm = x*normlize
    return (x_norm[0]**2+x_norm[1]/1000-11)**2+(x_norm[0]+(x_norm[1]/1000)**2-7)**2

normlize = np.array([1, 1000])
x_center = np.array([-1,-1])
step = 0.5
x0 = np.vstack((x_center, x_center+np.diag((step,step))))
xatol,fatol = 1e-3,1e-3
res = minimize(cost_function2,x_center,method='nelder-mead',options={'initial_simplex':x0,'xatol': xatol,'fatol': fatol, 'disp': True})
print(res.x)
x_ = normlize[0]*np.array([x[0] for x in xs])
y_ = normlize[1]*np.array([x[1] for x in xs])
ys = np.array([f2(x,y) for x,y in zip(x_,y_)])

fig, ax = plt.subplots()
c = ax.pcolormesh(xx,yy,z.T,cmap='jet')
fig.colorbar(c, ax=ax)
ax.plot(x_,y_,'r',x_[0],y_[0],'go',x_[-1],y_[-1],'y+',markersize=6)

fig2 = plt.figure()
ax1 = plt.subplot(221)
ax2 = plt.subplot(222)
ax3 = plt.subplot(212)
ax1.plot(x_)
ax1.set_title('x')
ax2.plot(y_)
ax2.set_title('y')
ax3.plot(ys)

结果：

Optimization terminated successfully.
         Current function value: 0.000002
         Iterations: 37
         Function evaluations: 69
[-3.77937271 -3.28299281]