分段函数拟合断点连续且一阶导数连续

Question

请问分段函数用python-lm拟合的时候，怎么添加条件让断点处连续且一阶导数连续？

Answer 1

在使用Python的 lmfit 模块进行分段函数拟合时，可以添加条件来强制断点处连续且一阶导数连续。具体操作如下：

1. 导入 lmfit 模块，并定义分段函数。

import numpy as np
from lmfit.models import ExpressionModel

def piecewise_linear(x, x0, y0, k1, k2):
    return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

该分段函数是由两个斜率分别为 k1 和 k2 的直线拼接而成，且在 x=x0 处分界。

1. 创建 lmfit.Model 对象并定义拟合参数。在拟合参数中，向模型参数中添加 min=x0, max=x0 来限制 x0 的取值范围，并定义 k2 为 k1 的别名。同时，定义两个约束函数，分别表示在 x=x0 处连续和一阶导数连续。

model = ExpressionModel('a*np.piecewise(x, [x < c], [lambda x:b*x + e-b*c, lambda x:d*x + e-d*c])',
            independent_vars=['x'])

params = model.make_params(a=1, b=0, c=50, d=0, e=0)
params.add('k1', expr='b')
params.add('k2', expr='k1')
params['c'].min = 0
params['c'].max = 100

def continuity(params):
    c = params['c'].value
    k1 = params['k1'].value
    k2 = params['k2'].value
    e = params['e'].value
    return [k1*c + e-k2*c]

def diff_continuity(params):
    c = params['c'].value
    k1 = params['k1'].value
    k2 = params['k2'].value
    return [k1-k2]
    
model.set_param_hint('k2', expr='k1')
model.set_param_hint(name='continuity', value=continuity)
model.set_param_hint(name='diff_continuity', value=diff_continuity)

其中，continuity 和 diff_continuity 分别为两个约束函数。continuity 函数返回一个列表，其中包含一个值，即在 x=x0 处断点处的连续条件公式的值。diff_continuity 函数返回一个列表，其中包含一个值，即在 x=x0 处一阶导数连续条件公式的值。这两个函数通过 lmfit 中的 model.set_param_hint 方法来定义。

1. 调用 lmfit.minimize 函数进行拟合。在调用该函数之前，将上面定义的约束函数加入到 fit_kws 参数中，代码如下：

result = model.fit(y, params, x=x, fit_kws={'constraints': [continuity, diff_continuity]})

fit_kws 参数是 lmfit.minimize 函数的关键字参数，通过该参数传递额外的拟合约束条件。

完整的拟合代码如下：

import numpy as np
from lmfit.models import ExpressionModel

# 分段函数定义
def piecewise_linear(x, x0, y0, k1, k2):
    return np.piecewise(x, [x < x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

# 创建 lmfit.Model 对象，定义拟合参数和约束条件函数
model = ExpressionModel('a*np.piecewise(x, [x < c], [lambda x:b*x + e-b*c, lambda x:d*x + e-d*c])',
            independent_vars=['x'])

params = model.make_params(a=1, b=0, c=50, d=0, e=0)
params.add('k1', expr='b')
params.add('k2', expr='k1')
params['c'].min = 0
params['c'].max = 100

def continuity(params):
    c = params['c'].value
    k1 = params['k1'].value
    k2 = params['k2'].value
    e = params['e'].value
    return [k1*c + e-k2*c]

def diff_continuity(params):
    c = params['c'].value
    k1 = params['k1'].value
    k2 = params['k2'].value
    return [k1-k2]
    
model.set_param_hint('k2', expr='k1')
model.set_param_hint(name='continuity', value=continuity)
model.set_param_hint(name='diff_continuity', value=diff_continuity)

# 调用 lmfit.minimize 进行拟合
result = model.fit(y, params, x=x, fit_kws={'constraints': [continuity, diff_continuity]})