AIDROID

Subtitle

0%

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

Posted on Edited on In , Views: Valine:

原文链接：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方法，函数的参数也会有相应的限制，这在具体的方法介绍中具体介绍。

1
2
3
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方法访问|

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

1
2
3
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单纯形法，该方法通过对可行域顶点不断迭代选出最优解。由于是迭代寻优策略，可以指定寻优范围。

2.2.2 powell

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

2.2.3 bfgs

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

2.2.4 newton-cg

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

2.2.5 trust-constr

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

三、实例求解

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两组初始假设值用于比较初始值对求解的影响。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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]])

        各个方法的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# 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的参数。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# 有限制的最优化
# 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)

​