优化 | 利用SciPy求解非线性规划问题

作者：莫斑炜

编者按：本文使用SciPy的optimize模块来求解非线性规划问题，结合实际例子，引入非线性规划问题的求解算法及相应函数的调用。

本文提纲

• 一维搜索／单变量优化问题

• 无约束多元优化问题

• 非线性最小二乘问题

• 约束优化问题

非线性规划问题的目标函数或约束条件是非线性的。本文使用SciPy的optimize模块来求解非线性规划问题。

目标函数和约束条件是否连续光滑是非常重要的性质，这是因为如果光滑，则所有决策变量可微，多变量函数的偏导数组成的向量为梯度，梯度是指向目标函数增长最快的方向。将目标函数梯度作为搜索方向，对非线性规划问题的求解具有重要的意义。这些函数或其导数\梯度的不连续性给许多现有的非线性优化问题的求解带来了困难。在下文中，我们假设这些函数是连续且光滑的。

# Importing Modules
from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np
import sympy

1、一维搜索／单变量优化问题（Univariate Optimization）

无约束非线性规划最简单的形式是一维搜索。一维搜索通常作为多维优化问题中的一部分出现，比如梯度下降法中每次最优迭代步长的估计。求解一维搜索常用的两类方法是函数逼近法和区间收缩法。其中函数逼近法是指用较简单的函数近似代替原来的函数，用近似函数的极小点来估计原函数的极小点，比如牛顿法；区间收缩法对于一个单谷函数通过迭代以不断缩小该区间的长度，当区间长度足够小时，可将该区间中的一点作为函数的极小点，比如黄金分割法。

e.g. 最小化一个单位体积的圆柱体的表面积。

Objective              

   s.t.                      

目标函数为圆柱体的表面积，约束条件为圆柱体体积为1。变量为r、h，这是一个带等式约束的二维优化问题。根据等式约束条件，可得              ，代入目标函数得到一维搜索：              。该问题可通过求解导数为0的方法来求解。

r, h = sympy.symbols("r, h")
Area = 2 * sympy.pi * r**2 + 2 * sympy.pi * r * h
Volume = sympy.pi * r**2 * h
h_r = sympy.solve(Volume - 1)[0]
Area_r = Area.subs(h_r)
rsol = sympy.solve(Area_r.diff(r))[0]
rsol.evalf()
# 0.541926070139289
 
 
# 再验证二阶导数为正且rsol对应最小值
Area_r.diff(r, 2).subs(r, rsol)
Area_r.subs(r, rsol)
_.evalf()
# 5.53581044593209

使用SciPy中optimize模块中的brent()函数来最小化一维函数。brent方法是一种混合方法，结合了牛顿法和二分法，既能保证稳定性又能快速收敛，通常是SciPy一维搜索的首选方法。

首先定义一个Python函数f作为目标函数。将f传给optimize.brent，参数brack表示指定算法的开始区间。

def f(r):
    return 2 * np.pi * r**2 + 2 / r
r_min = optimize.brent(f, brack=(0.1, 4))
r_min
# 0.5419260772557135
f(r_min)
# 5.535810445932086

通过以上求解结果可看出，r取值约为0.54，目标函数取值约为5.54。注意进行优化之前，最好通过目标函数可视化以确定合适的搜索区间或起点。

2、无约束多元优化问题

无约束多元优化问题的求解要比一维搜索问题的求解困难得多。多变量情况下，求解非线性梯度方程根的解析方法是不可行的，而黄金分割搜索法中使用的区间形式也不能直接用。最基本的方法是考虑目标函数f(x)在给定点x处的梯度。负梯度方向是指向函数f(x)减少最多的方向，即负梯度方向是函数的最速下降方向。牛顿多元优化方法是对最速下降法的一种改进，它可以提高收敛性。在单变量情况（一维搜索）下，牛顿法可以看作函数的局部二次逼近。在多变量情况下，与最速下降法相比，步长被函数Hessian矩阵的逆的所代替。

SciPy库中，函数optimize.fmin_ncg使用牛顿法，optimize.fmin_ncg需要的参数包括目标函数和搜索起点，需要用于计算梯度的函数及用于计算Hessian矩阵的函数。

e.g.              

使用牛顿法进行求解，需要知道梯度和Hessian矩阵，对每个变量使用sympy.diff函数（用于求导）以获得梯度和Hessian矩阵。

x1, x2 = sympy.symbols("x_1, x_2")
f_sym = (x1-1)**4 + 5 * (x2-1)**2- 2*x1*x2
fprime_sym = [f_sym.diff(x_) for x_ in (x1, x2)]
# Gradient
sympy.Matrix(fprime_sym)
fhess_sym = [[f_sym.diff(x1_, x2_) for x1_ in (x1, x2)] for x2_ in (x1, x2)]
# Hessian
sympy.Matrix(fhess_sym)
 
 
# 构造函数表达式
f_lmbda = sympy.lambdify((x1, x2), f_sym, 'numpy')
fprime_lmbda = sympy.lambdify((x1, x2), fprime_sym, 'numpy')
fhess_lmbda = sympy.lambdify((x1, x2), fhess_sym, 'numpy')
 
 
def func_XY_to_X_Y(f):
    """
    Wrapper for f(X) -> f(X[0], X[1])
    """
    return lambda X: np.array(f(X[0], X[1]))
f = func_XY_to_X_Y(f_lmbda)
fprime = func_XY_to_X_Y(fprime_lmbda)
fhess = func_XY_to_X_Y(fhess_lmbda)
# f、fprime、fhess不仅作用于单个值，还作用于整个向量
 
 
# Newton method
x_opt = optimize.fmin_ncg(f, (0, 0), fprime=fprime, fhess=fhess) # (0,0)为搜索起点
x_opt
# array([1.88292613, 1.37658523])

       

目标函数极小值用红星标出

由于牛顿法需要提供梯度和Hessian矩阵的计算函数，有时它们很难计算。下面介绍的方法可以不用传入计算梯度和Hessian矩阵的函数。用optimize.fmin_bfgs（BFGS算法）和optimize.fmin_cg（共轭梯度法）求解，不必为Hessian函数提供函数。如果没有计算梯度的函数也可以不传入，直接调用。

# 仅提供梯度，不提供hessian矩阵
# BFGS method
x_opt = optimize.fmin_bfgs(f, (0, 0), fprime=fprime)
x_opt
 
 
# conjugate gradient method
x_opt = optimize.fmin_cg(f, (0, 0), fprime=fprime)
x_opt

也可以在不提供梯度函数的情况下使用：

x_opt = optimize.fmin_bfgs(f, (0, 0))

但是如果传递了计算梯度的函数，求解结果将更好，求解速度也更快。

在梯度和Hessian矩阵都未知的情况下，可使用BFGS算法（拟牛顿方法）。如果梯度和Hessian矩阵都是已知的，那么牛顿法是通常收敛速度最快的方法，但是对于大规模问题牛顿法的每次迭代的计算量较大。虽然BFGS算法和共轭梯度法在理论上的收敛速度比牛顿法慢，但它们有时稳定性更好。

以上介绍的多元优化方法只能收敛到局部最优，即使存在全局最优解，这些方法也很难跳出局部最优，对于这种情况可以通过蛮力搜索找到更好的搜索起点来提升算法的性能。SciPy提供optimize.brute函数执行这种系统搜索。

e.g.              

def f(X):
    x, y = X
    return (4 * np.sin(np.pi * x) + 6 * np.sin(np.pi * y)) +  (x- 1)**2 + (y- 1)**2

调用optimize.brute，目标函数f作为第一个参数，切片对象元组作为第二个参数，每个坐标一个。切片对象指定了搜索范围，并通过设置参数finish=none直接给出蛮力的搜索方案。

x_start = optimize.brute(f, (slice(-3, 5, 0.5),  slice(-3, 5, 0.5)), finish=None)
x_start
# array([2. , 1.5])
f(x_start)
# 得到更好的搜索起点后代入optimize.fmin_bfgs中
x_opt = optimize.fmin_bfgs(f, x_start)
# 可以和任意搜索起点作对比
x_opt = optimize.fmin_bfgs(f,(0,0))

可以看出目标函数更优，且求解时间更短，迭代次数更少，速度更快。

minimize()函数为scipy.optimize中的多变量标量函数提供了无约束和约束最小化算法的通用接口。SciPy提供为多变量优化求解器提供统一的接口optimize.minimize()，它通过method参数为求解器指定特定的函数，这样可以更容易地在不同的求解器之间切换。一维搜索的接口为optimize.scalar_minimize。

x_opt = optimize.fmin_bfgs(f, x_start)  # 等价于下面的语句
result = optimize.minimize(f, x_start, method= 'BFGS') # 也可不指定method
x_opt = result.x

3、非线性最小二乘问题

                   

optimize.leatsq函数用于求解最佳最小二乘拟合。

e.g.              

# 观测值
beta = (0.25, 0.75, 0.5)
def f(x, b0, b1, b2):
    return b0 + b1 * np.exp(-b2 * x**2)
xdata = np.linspace(0, 5, 50)
y = f(xdata, *beta)
# 加入噪音
ydata = y + 0.05 * np.random.randn(len(xdata))
 
 
def g(beta):
    return ydata- f(xdata, *beta)
beta_start = (1, 1, 1)
beta_opt, beta_cov = optimize.leastsq(g, beta_start)
beta_opt
# array([0.24935676, 0.74672532, 0.49918151])
# 求解结果接近beta(0.25, 0.75, 0.5)

拟合曲线：

             

SciPy的optimize模块还通过函数optimize.curve_fit为非线性最小二乘拟合提供了另一种接口，它不用向函数传递一个带剩余误差(residuals)的显式函数。

beta_opt, beta_cov = optimize.curve_fit(f, xdata, ydata)
beta_opt
# array([0.24935676, 0.74672532, 0.49918151])

4、约束优化问题

约束条件增加了非线性规划问题求解的复杂性。可以使用拉格朗日乘子法（Lagrange multiplier），通过引入附加变量将约束优化问题转化为无约束问题。

e.g.

             

求构成的最大体积。

求解该问题使用拉格朗日乘子法，拉格朗日函数为              。可以找到使得该函数梯度为0的驻点求解该问题。

x = x0, x1, x2, l = sympy.symbols("x_0, x_1, x_2, lambda")
f = x0 * x1 * x2
g = 2 * (x0 * x1 + x1 * x2 + x2 * x0)- 1
L = f + l * g
grad_L = [sympy.diff(L, x_) for x_ in x]
sols = sympy.solve(grad_L)
sols
# 求解结果中负值舍去
g.subs(sols[0])
# 0
f.subs(sols[0])

SciPy的optimize模块中的函数slsqp()或通过optimize.minimize()函数设置method为'SLSQP'，为了用SciPy的slsqp求解器求解问题，我们需要为目标函数和约束函数定义python函数。由于optimize模块优化函数用于解决最小化问题，而这里是求解最大化问题，故将函数f设置为原始目标函数的相反数。

def f(X):
    return -X[0] * X[1] * X[2]
def g(X):
    return 2 * (X[0]*X[1] + X[1] * X[2] + X[2] * X[0])- 1

接下来定义g(x)=0的约束。通过字典形式传入约束条件，该字典中key（value）是type（'eq'或'ineq'）、fun（约束条件）、jac（约束函数的Jacobian矩阵）和args（约束函数的其他参数和计算其Jacobian的函数）。最后调用optimize.minimize函数。

constraint = dict(type='eq', fun=g)
result = optimize.minimize(f, [0.5,1,1.5],method='SLSQP',constraints=[constraint])
result

可以看出通过以上方法求得的解与通过使用拉格朗日乘子法得到的解基本一致。

要解决不等式约束问题，只需在约束字典中设置type='ineq'，并给出相应的不等式函数。

e.g.

             

def f(X):
    return (X[0]- 1)**2 + (X[1]- 1)**2
def g(X):
    return X[1]- 1.75- (X[0]- 0.75)**4
constraints = [dict(type='ineq', fun=g)]
x_opt = optimize.minimize(f, (0, 0), method='BFGS').x  #无约束
x_cons_opt = optimize.minimize(f(0,0),method='SLSQP',constraints=constraints).x

       

约束问题的可行域为灰色阴影，红星和蓝星分别为有约束和无约束问题的最优解。

参考文献：

Johansson R. Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib[M]. Apress, 2019.

备注：公众号菜单包含了整理了一本AI小抄，非常适合在通勤路上用学习。

往期精彩回顾
 
 
 
 
2019年公众号文章精选适合初学者入门人工智能的路线及资料下载机器学习在线手册深度学习在线手册AI基础下载（第一部分）备注：加入本站微信群或者qq群，请回复“加群”加入知识星球（4500+用户，ID：92416895），请回复“知识星球”

喜欢文章，点个在看