python模块：Scipy.optimize.minimize规划问题求解

一、模块介绍

1.1模块功能

Scipy.optimize是Scipy中一个用于解决数学模型中优化类模型的子包，该子包中又包含了多个子功能模块见下表，不同方法不同条件求解最优化模型。本节介绍minimize对一般规划问题的模型建立与求解。

问题类型	模块
多元标量函数的有/无约束最小化	minimize
最小二乘法最小化	least_squares
单变量函数最小化器	minimize_scalar
线性规划	linprog

1.2模型介绍

多元标量函数的最小化，是数学规划模型中更为一般的模型，该模块包括有限制性约束和无限制性约束的最小化，而对于限制性约束又分为线性约束和非线性约束。这种更为一般的模型需要针对具体的问题假设选择特定的方法进行求解。

在数学规划模型中，minimize提供的方法能够解决无/有（线性、非线性）约束的多个决策变量目标函数的最优化问题，但是由于该模块是依据函数导数与梯度进行求解，不能够求解整数规划、01规划等问题。

二、模块源分析与参数解释

2.1模块整体介绍

对于整个多元标量的最小化问题，主要的api为optimize包下的minimize函数。该minimize函数的完整参数详情如下。同时，对于不同的method方法，函数的参数也会有相应的限制，这在具体的方法介绍中具体介绍。

minimize(fun, x0, args=(), method=None, jac=None, hess=None,
             hessp=None, bounds=None, constraints=(), tol=None,
             callback=None, options=None)

该函数主要参数意义解释：

fun	待求解的目标函数
x0	初始猜测的数组
args	目标函数带参数时需要指定
method	选择的方法
jac	jacobian矩阵或梯度函数（梯度求解的方法需要传入）
hess	hessian矩阵
hessp	黑塞向量积
bounds	寻优范围
constraints	限制约束
tol	浮动，可以理解为单次寻找的步长上限，越小精度越高
callback	回调函数
options	选项字典，不同方法有不同选项

完整的方法options可通过optimize的show_options方法访问

from scipy.optimize import show_options
print(show_options(method=’minimize’)

该模块提供求解方法字典如下, 通过optimize包导入minimize后，对minimize传入method参数指定方法。而对于不同方法的选用,便有不同的求解策略。我们选择部分主要方法介绍。

MINIMIZE_METHODS = ['nelder-mead', 'powell', 'cg', 'bfgs', 'newton-cg',
                    'l-bfgs-b', 'tnc', 'cobyla', 'slsqp', 'trust-constr',
                    'dogleg', 'trust-ncg', 'trust-exact', 'trust-krylov']

无约束多元标量最小化选择以下4个方法：

nelder-mead
powell
bfgs
newton-cg

有约束多元标量最小化选择方法：

trust-constr

对于更多的方法可以访问该文档查询：完整方法文档

2.2 各方法详细介绍

2.2.1 nelder-mead

方法完名称为Nelder-Mead Simplex algorithm Nelder-Mead单纯形法，该方法通过对可行域顶点不断迭代选出最优解。由于是迭代寻优策略，可以指定寻优范围。

该方法对应参数：
scipy.optimize.minimize(fun, x0, args=(), method=’Nelder-Mead’, bounds=None, tol=None, callback=None, options={‘func’: None, ‘maxiter’: None, ‘maxfev’: None, ‘disp’: False, ‘return_all’: False, ‘initial_simplex’: None, ‘xatol’: 1e-5, ‘fatol’: 1e-5, ‘adaptive’: False})

特殊参数：
bounds 该方法不要求函数连续，可指定寻优范围
特殊选项：
xatol/fatol 解与函数值迭代之间绝对误差上限值，值越小求解越慢，精度越高

2.2.2 powell

方法名称为鲍威尔算法，该方法通过共轭方向加速收敛，其它与单纯形法类似。

该方法对应参数：
scipy.optimize.minimize(fun, x0, args=(), method=’Powell’, bounds=None, tol=None, callback=None, options={‘func’: None, ‘xtol’: 0.0001, ‘ftol’: 0.0001, ‘maxiter’: None, ‘maxfev’: None, ‘disp’: False, ‘direc’: None, ‘return_all’: False})

特殊参数：
bounds 该方法不要求函数连续，可指定寻优范围
特殊选项：
xatol/fatol 解与函数值迭代之间绝对误差上限值，值越小求解越慢，精度越高

2.2.3 bfgs

方法名称为 BFGS（Broyden-Fletcher-Goldfarb-Shanno algorithm）拟牛顿法。算法利用梯度下降进行收敛，因此不可指定范围。该方法未传入梯度函数时通过第一次的差值估计。

该方法对应参数：
scipy.optimize.minimize(fun, x0, args=(), method=’BFGS’, jac=None, tol=None, callback=None, options={‘gtol’: 1e-05, ‘norm’: inf, ‘eps’: 1.4901161193847656e-08, ‘maxiter’: None, ‘disp’: False, ‘return_all’: False, ‘finite_diff_rel_step’: None})

特殊参数：
         jac        梯度函数
特殊选项：
         gtol       梯度范数上限
norm      范数浮动上限
eps 求解浮动上限

2.2.4 newton-cg

方法名称为NCG （Newton-Conjugate-Gradient algorithm）牛顿共轭梯度算法，该方法通过梯度函数，二阶导数的hessian矩阵形式将目标函数拟合为二次函数形式，通过该二次式梯度下降求解。

该方法对应参数：
scipy.optimize.minimize(fun, x0, args=(), method=’Newton-CG’, jac=None, hess=None, hessp=None, tol=None, callback=None, options={‘xtol’: 1e-05, ‘eps’: 1. e-08, ‘maxiter’: None, ‘disp’: False, ‘return_all’: False})

特殊参数：
         jac        梯度函数
         hess         hessian矩阵
         hessp hessian向量积（与hess只需要定义一个）
特殊选项：
         xtol        解的平均相对误差上限
eps 求解浮动上限

2.2.5 trust-constr

方法名称为Trust-Region Constrained Algorithm信赖域约束算法，该方法需要用用约束模块定义线性和非线性约束。线性约束与linprog中约束定义类似，而非线性约束需要定义jacobian矩阵和hessian矩阵。

该方法对应参数：
scipy.optimize.minimize(fun, x0, args=(), method=’trust-constr’, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options={‘grad’: None, ‘xtol’: 1e-08, ‘gtol’: 1e-08, ‘barrier_tol’: 1e-08, ‘sparse_jacobian’: None, ‘maxiter’: 1000, ‘verbose’: 0, ‘finite_diff_rel_step’: None, ‘initial_constr_penalty’: 1.0, ‘initial_tr_radius’: 1.0, ‘initial_barrier_parameter’: 0.1, ‘initial_barrier_tolerance’: 0.1, ‘factorization_method’: None, ‘disp’: False})

特殊参数：
        bounds 变量定义域
         jac        梯度函数
         hess         hessian矩阵
         hessp hessian向量积（与hess只需要定义一个）
constraints 线性约束和非线性约束组成的约束列表
特殊选项：
        xatol/fatol   解与函数值迭代之间绝对误差上限值，值越小求解越慢，精度越高

三、实例求解

3.1 实例模型与知识补充

3.1.1Rosenbrock函数

rosenbrock是一个常用于检测优化算法效果的函数Rosenbrock函数_百度百科 (baidu.com)。

该函数可以用以形式表示：

$f(X)=\sum_{i=1}^{N-1}100(x_{i+1}-x_{i}^{2} )^{2}+(1-x_{i})^{2}$

其中，X为N为向量表示多元标量

$X=(x_{1},...,x_{N})$

该函数的全局最优解为： $X_{g}=(1,1,...,1)$

3.1.2 jacobian矩阵和hessian矩阵

这两个矩阵是用于梯度求解方法，jacobian矩阵是函数一阶偏导形成的一个向量。hessian矩阵是二阶偏导后得到的一个N阶方阵，当指定X=X0时，在该点的hessian矩阵可表示为：

通过这两个矩阵可以将目标函数通过二次函数形式拟合。

对于数学部分了解以上部分即可，具体的数学解释可以参考：黑塞矩阵_百度百科 (baidu.com)

3.2 实例求解

3.2.1无约束多元标量最小化

根据上面的介绍，我们以rosenbrock函数作为实例演示（下称rosen函数）。

先定义rosen的目标函数，梯度函数和hessian矩阵以及一些假设变量如假定范围为bounds。设定x1，x2两组初始假设值用于比较初始值对求解的影响。

from scipy.optimize import minimize
import numpy as np

# 两组初始迭代值
x1=np.array([1,1.2,0.9,1.3,0.5])
x2=np.array([1,1.2,0.9,1.3,0.3])
# 以rosenbrock函数来评价优化方法
def rosen(x):
    return sum(100*(x[1:]-x[:-1]**2.0)**2+(1-x[:-1])**2.0)
# 带参的rosen函数
def rosen_(x,a,b):
    return sum(a*(x[1:]-x[:-1]**2.0)**2+(1-x[:-1])**2.0)+b
# rosen函数的梯度
def rosen_der(x):
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
    der[-1] = 200*(x[-1]-x[-2]**2)
    return der
# rosen函数的hessian矩阵
def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H
# 限定寻优范围
bounds=np.array([[0,2],[0,2],[0,2],[0,2],[None,None]])

各个方法的实现如下：

# nelder-mead单纯形方法 初始值依赖
res11=minimize(rosen_,x1,args=(100,1),method='nelder-mead',options={'xatol': 1e-8, 'disp':True})
res12=minimize(rosen_,x1,args=(100,1),method='nelder-mead',options={'xatol': 1e-8, 'disp':True})
print('========================')
print(res11)
print('========================')
print(res12)


# powell方法
res21=minimize(rosen_,x1,args=(100,1),method='powell',tol=1e-5,bounds=bounds,options={'disp':True})
res22=minimize(rosen_,x2,args=(100,1),method='powell',bounds=bounds,options={'disp':True})
# ftol和xtol分别为函数值与x解的相对误差上限
res21_=minimize(rosen_,x1,args=(100,1),method='powell',options={'ftol':1e-8,'xtol':1e-8,'disp':True})
print('初始迭代值会对算法的收敛速度和精度产生影响========================')
print(res21.x)
print(res22.x)
print('对于精度过低，可以指定浮动上限值，通过增加代数提高算法精度========================')
print(res21_.x)


# bfgs拟牛顿方法  由于拟牛顿法采用梯度计算，不能够指定变量的寻优范围 gtol为梯度范数上限
res31=minimize(rosen,x1,method='bfgs',options={'gtol':1e-8,'disp':True})
res32=minimize(rosen,x1,method='bfgs',jac=rosen_der,options={'gtol':1e-8,'disp':True})
print('jac参数为函数梯度参数，若不指定，则通过差值估计,两者效果相差不大============')
print(res31)
print(res32)

# ncg 牛顿共轭梯度方法，先求解hessian矩阵得到二次函数拟合目标函数，再通过梯度下降求解最优
res41=minimize(rosen,x1,method='Newton-cg',jac=rosen_der,hess=rosen_hess,options={'xtol':1e-8,'disp':True})
print(res41)

3.2.2 有约束的多元标量最小化

我们选择二元标量的rosen函数，并添加约束，得到一个优化模型：

$min 100(x_{1}-x_{0}^{2})^{2}+(1-x_{0})^{2}$

s.t.

$\begin{cases} & \text{ } x_0+2x_1\leq 1 \\ & \text{ } x_0^2+x_1\leq 1 \\ & \text{ } x_0^2-x_1\leq 1 \\ & \text{ } 2x_0+x_1 = 1 \\ & \text{ } 0 \leq x_0 \leq 1 \\ & \text{ } -0.5 \leq x_1 \leq 2 \end{cases}$

按照下步骤建立模型并求解，其中线性约束部分与linprog模块类似，通过矩阵表示多个线性约束条件。而非线性约束需要定义对应的约束求解函数cons_f，jacobian矩阵函数和hessian矩阵的线性组合函数cons_J/cons_H.如下所示。两个约束需要通过LinearConstraint和NonlinearConstraint封装为对应的结构，并组合成列表形式作为constraints的参数。

# 有限制的最优化
# Trust-Region Constrained Algorithm信赖域约束算法
# 定义一个带有线性和非线性约束的优化函数
from scipy.optimize import Bounds,LinearConstraint,NonlinearConstraint
# bounds为定义域，参数分别为变量的下限和上限向量表示形式
bounds_t=Bounds([0, -0.5], [1.0, 2.0])
# 与linprog中定义线性约束相同，通过矩阵和向量表示参数
linear_constraint=LinearConstraint([[1,2],[2,1]],[-np.inf,1],[1,1])
# 非线性约束需要jacobian矩阵和hessian矩阵辅助定义
# 返回非线性约束原函数
def cons_f(x):
    return [x[0]**2 + x[1], x[0]**2 - x[1]]
# 对应返回jacboian矩阵
def cons_J(x):
    return [[2*x[0], 1], [2*x[0], -1]]
# 返回hessian矩阵的线性组合形式
def cons_H(x,v):
    return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]])
from scipy.optimize import NonlinearConstraint
nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
# 约束条件
x3=np.array([0.5,0])
res5=minimize(rosen, x3, method='trust-constr', jac=rosen_der,constraints=[linear_constraint, nonlinear_constraint],hess=rosen_hess,options={'verbose': 1}, bounds=bounds_t)
print(res5)

四、参考

[1] Scipy v1.9.0版用户指南 Optimization (scipy.optimize) — SciPy v1.9.0 Manual

[2] osego 指南优化与寻根 (scipy.optimize ) — SciPy v1.8.0.dev0+1869.838cfbe Manual (osgeo.cn)

文章出处登录后可见！

已经登录？立即刷新