scipy.optimize应用-使用curve_fit和minimize进行函数拟合与非线性规划 Python

艽野尘梦better

已于 2023-07-11 13:52:12 修改

阅读量4.1k

点赞数 5

分类专栏： Python应用文章标签： scipy python 开发语言

于 2023-07-11 13:48:56 首次发布

本文链接：https://blog.csdn.net/qq_45270849/article/details/131416108

版权

Python应用专栏收录该内容

8 篇文章 1 订阅

订阅专栏

3.3 基于拟合数据得到的函数求解非线性规划

引言

接上一篇文章，我们使用自有数据拟合曲面表达式并绘制3d曲面图：如何根据自有数据拟合绘制3D曲面图-Python matplotlib_python 拟合三维曲线_艽野尘梦better的博客-CSDN博客https://blog.csdn.net/qq_45270849/article/details/130945731?spm=1001.2014.3001.5502

在该文中，我们使用scipy.optimize.curve_fit来拟合自有数据的函数，获得函数中未知参数的值。在本文中，我们将结合scipy.optimize.curve_fit和scipy.optimize.minimize两个方法实现拟合自有数据函数表达式，并进行非线性规划寻找拟合后函数在带约束条件下的最大值或最小值。

主要参考如下：

scipy.optimize.minimize — SciPy v1.11.1 Manualhttps://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html 非线性规划（scipy.optimize.minimize） - 简书 (jianshu.com)https://www.jianshu.com/p/94817f7cc89b Python小白的数学建模课-12.非线性规划 - 知乎 (zhihu.com)https://zhuanlan.zhihu.com/p/393509520

一、数据描述

为了方便演示，我们生成两个自变量x1和x2，一个带噪声的因变量y，他们在理想上的函数形式是 $y=x_{1}^{2}-2x_{1}x_{2}+x_{2}^{2}$ ，由于y带噪声，所以我们要拟合的函数形式是 $y=ax_{1}^{2}+bx_{1}x_{2}+cx_{2}^{2}$ ，其中a、b、c是未知参数。为了贴合一般情形，我将数据放在一个dataframe中。

import numpy as np
import pandas as pd
from scipy.optimize import curve_fit
from scipy.optimize import minimize

x1=np.random.randint(-50,50,100)
x2=np.random.randint(-50,50,100)
y=x1**2-2*x1*x2+x2**2+np.random.normal(loc=0,scale=5,size=100)
data=pd.DataFrame([x1,x2,y],index=['x1','x2','y']).T
data.head()

二、拟合函数

2.1 curve_fit函数定义

为了方便理解，先给出curve_fit方法的定义：

scipy.optimize.curve_fit(f, xdata, ydata, p0=None, sigma=None, absolute_sigma=False, check_finite=None, bounds=(-inf, inf), method=None, jac=None, *, full_output=False, nan_policy=None, **kwargs)

其中：

f是要拟合的函数形式，函数的参数是自变量和需要确定的参数；

xdata是自变量数据，根据自变量个数，如果是一个则为一维数组，如果多于一个，则为二维数组，数组行数是自变量个数；

ydata是因变量数据，应该为一维数组，长度与xdata相同

p0是需要拟合的未知参数的可能取值

bounds是未知参数的范围，如果元组中是整数如(-5, 5)，则所有未知参数均在此范围内寻找，也可以分别制定每个未知参数的范围，形式如([-5,-5,-5],[5,5,5])

method是用于优化的方法，有{‘lm’, ‘trf’, ‘dogbox’}三种取值，对于没有约束的问题使用‘lm’，对于有约束即指定了bounds的问题使用‘trf’，如果数据量小于未知参数个数，则‘lm’将不起作用，此时使用‘trf’，‘dogbox’。

return：主要使用popt和pcov两个返回值，其中popt是使误差平方和最小化的最优拟合参数值，pcov是popt的估计近似协方差，更多参数和返回值定义详见函数的官方文档。

2.2 确定未知参数

基于以上基础，我们可以拟合出函数中未知参数a、b、c的值：

def func(X,a,b,c):
    y=a*X[0]**2+b*X[0]*X[1]+c*X[1]**2
    return y
x=np.array(data[['x1','x2']].T).reshape(2,-1)
y=np.array(data['y'])
popt, pcov =curve_fit(func,x,y)
popt

三、函数寻优

在scipy.optimize中提供了一系列用于函数寻优的方法，其中scipy.optimize.brent() 用于求解单变量无约束优化问题；scipy.optimize.fmin() 用于求解多变量无约束优化问题；求解非线性规划问题则使用scipy.optimize.minimize()。本文使用minimize求解非线性规划，我们假设寻优问题如下，其中a、b、c是之前拟合的未知参数：

$min\ ax_1^{2}+bx_1x_2+cx_2^{2}\\s.t.\begin{cases} x_1\ge 2x_2\\ 1\le x_1\le 5\\ 1\le x_2\le 5\\ \end{cases} \\ \\$

3.1 minimize函数定义

同样先给出minimize的函数定义：

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

其中：

fun: callable f(x,*args) 是目标函数，以函数形式表示，可以通过 *args 传递参数。

x0: nadarray, shape(n,) 搜索算法的初始值，n是决策变量个数。必须给出

args: tuple 可选项，将可变参数传递给目标函数 fun、导函数 jac和二阶导函数hess。

method: str或callable 可选项，选择优化算法。默认算法为 BFGS, L-BFGS-B, SLSQP（取决于问题有没有边界条件和约束条件）

jac: 可选项，梯度计算方法。可选择{callable, ‘2-point’, ‘3-point’, ‘cs’, bool}。该选项只能用于 CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact 和 trust-constr 算法。

hess: 可选项，Hessian 矩阵计算方法。可选择{callable, ‘2-point’, ‘3-point’, ‘cs’, HessianUpdateStrategy}。该选项只能用于 Newton-CG, dogleg, trust-ncg, trust-krylov, trust-exact 和 trust-constr 算法。

bounds: 可选项，变量的边界条件（上下限，lb<=x<=ub）。该选项只能用于 Nelder-Mead, L-BFGS-B, TNC, SLSQP, Powell 和 trust-constr 算法。

constraints: 可选项，定义约束条件如 f(x)>=0。该选项只能用于 COBYLA, SLSQP 和 trust-constr 算法，注意不同算法中对于约束条件的定义是不同的。

return：OptimizeResult对象。常用该对象的x和fun属性，分别表示最优解和最优函数值。

3.2 minimize函数可用优化算法

一般求解非线性规划问题使用‘SLSQP’算法。其他优化算法的定义摘自知乎youcans的文章。

无约束问题优化算法

method=‘CG’ ：非线性共轭梯度算法，只能处理无约束优化问题，需要使用一阶导数函数。
method=‘BFGS’ ：BFGS 拟牛顿法，只能处理无约束优化问题，需要使用一阶导数函数。BFGS 算法性能良好，是无约束优化问题的默认算法。
method=‘Newton-CG’ ：截断牛顿法，只能处理无约束优化问题，需要使用一阶导数函数，适合处理大规模问题。
method=‘dogleg’ ：dog-leg 信赖域算法，需要使用梯度和 Hessian（必须正定），只能处理无约束优化问题，
method=‘trust-ncg’ ：采用牛顿共轭梯度信赖域算法，需要使用梯度和 Hessian（必须正定），只能处理无约束优化问题，适合大规模问题。
method=‘trust-exact’：求解无约束极小化问题的信赖域方法，需要梯度和Hessian（不需要正定）。
method=‘trust-krylov’：使用Newton-GLTR 信赖域算法度，需要使用梯度和 Hessian（必须正定），只能处理无约束优化问题，适合中大规模问题。

边界约束条件问题优化算法

method=‘Nelder-Mead’：下山单纯性法，可以处理边界约束条件（决策变量的上下限），只使用目标函数，不使用导数函数、二阶导数，鲁棒性强。
method=‘L-BFGS-B’ ：改进的 BFGS 拟牛顿法，L- 指有限内存，-B 指边界约束，可以处理边界约束条件，需要使用一阶导数函数。L-BFGS_B 算法性能良好，消耗内存量很小，适合处理大规模问题，是边界约束优化问题的默认算法。
method=‘Powell’：改进的共轭方向法，可以处理边界约束条件（决策变量的上下限）。
method=‘TNC’ ：截断牛顿法，可以处理边界约束条件

带有约束条件问题优化算法

method=‘COBYLA’ ：线性近似约束优化方法，通过对目标函数和约束条件的线性逼近处理非线性问题。只使用目标函数，不需要导数或二阶导数值，可以处理约束条件。
method=‘SLSQP’ ：序贯最小二乘规划算法，可以处理边界约束、等式约束和不等式约束条件。SLSQP 算法性能良好，是带有约束条件优化问题的默认算法。
method=‘trust-constr’ ：信赖域算法，通用的约束最优化方法，适合处理大规模问题。

3.3 基于拟合数据得到的函数求解非线性规划

基于以上定义，我们的非线性规划问题可以如下求解：

#取值范围
bound1=(1,5)
bound2=(1,5)
bounds=(bound1,bound2)
a,b,c=popt
#约束条件
def constraint1(x):  # 不等式约束
    return x[0]-2*x[1]
con = {'type': 'ineq', 'fun': constraint1}#type为ineq表示不等式约束，f(x)>=0；type为eq表示等式约束
cons = ([con]) 
#初始值
x0=np.array([5,5])
res = minimize(func, x0, args=(a,b,c),method='SLSQP', bounds=bounds, constraints=cons)
print(res.x,res.fun)