首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >当导数的值受到约束时,如何拟合多项式?

问当导数的值受到约束时,如何拟合多项式?
EN

Stack Overflow用户
提问于 2021-06-15 18:01:30
回答 1查看 60关注 0票数 1

对于x的某个取值范围,是否存在对多项式函数f(x)的导数的最大值有约束的解的简单方法

与回答this question一样,来自scipy.optimize的curve_fit可以处理单个系数上的约束,如下例所示:

代码语言:javascript
复制
def func(x, a, b, c, d):
  return a + b * x + c * x ** 2 + d * x ** 3

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

popt_cons, _ = curve_fit(func, x, y, bounds=([-np.inf, 2, -np.inf, -np.inf], [np.inf, 2.001, np.inf, np.inf]))
print(popt_cons)
>>> [-0.14331349  2.         -0.95913556  0.10494372]

但是,如果我想要最佳拟合多项式,其中存在对某个x值范围的加速度(二阶导数)的最大值的约束,该怎么办?这意味着,通过对函数进行两次集成,假设在0到10之间使用x,则对2*c + 6*d*x的值有一个约束。有没有方法可以做到这一点,或者我必须从头开始构建?

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2021-06-15 20:34:57

curve_fit方法不支持附加约束。但是,您可以实现非线性最小二乘问题。

代码语言:javascript
复制
min ||f(x, coeffs) - y||^2

s.t.    lb <= coeffs <= ub
     f''(x, coeffs) <= max_val for all x_lb <= x <= x_ub

并使用minimize解决它。下面是一个如何通过np.polyval和np.polyder完成此操作的示例

代码语言:javascript
复制
import numpy as np
from scipy.optimize import minimize

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

def objective(coeffs):
    return np.linalg.norm(np.polyval(coeffs, x) - y)

def constraint(coeffs, deriv_order, x_lb, x_ub, max_val):
    deriv_coeffs = np.polyder(coeffs, deriv_order)
    # Evaluate the derivative for all x_lb <= x <= x_ub
    deriv_value = np.polyval(deriv_coeffs, x[(x >= x_lb) & (x <= x_ub)])
    return -1.0*deriv_value + max_val

# Each inequality constraint has the form fun(x) >= 0
cons = [{'type': 'ineq', 'fun': lambda coeffs: constraint(coeffs, 2, 0.0, 3.0, 20)}]
bnds = [(-np.inf, np.inf), (2, 5.001), (-np.inf, np.inf), (-np.inf, np.inf)]

poly_degree = 3
res = minimize(objective, x0=2.0*np.ones(poly_degree+1), bounds=bnds, constraints=cons)

注意,每个不等式约束都有fun(x) >= 0的形式,即我们有-f''(x, coeffs) + max_val >= 0,我们使用x_lb = 0.0,x_ub = 3.0和max_val = 20作为二阶导数。最后,res.x包含多项式系数。

票数 2
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/67984130

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档