SciPy.optimize.least_squares() Objective 函数问题

Question

我试图通过优化三个未知参数 a、b 和 c0 来最小化高度非线性函数。我试图在 Python 3.

中复制赌场轮盘球的一些控制方程

这是研究论文的 link： http://www.dewtronics.com/tutorials/roulette/documents/Roulette_Physik.pdf

我将在论文中引用等式 (35) 和 (40)。基本上，我用秒表测量轮盘上旋转的轮盘球的圈数。对于每个连续的圈，由于非保守摩擦力的动量损失，圈时间将增加。然后，我使用等式 (40) 中的 Levenberg-Marquardt 最小二乘法进行这些时间测量并拟合等式 (35)。

我的问题有两个： (1)我用的是scipy.optimize.least_squares()方法='lm'，不知道objective函数怎么写！现在我已经完全按照论文中的方式编写了函数：

def fall_time(k,a,b,c0):

    F = (1 / (a * b)) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))

    return F

def parameter_estimation_function(x0,tk):

    a = x0[0]
    b = x0[1]
    c0 = x0[2]

    S = 0

    for i,t in enumerate(tk):

        k = i + 1

        S += (t - fall_time(k,a,b,c0))**2

    return [S,1,1]

sol = least_squares(parameter_estimation_function,[0.1,0.8,-0.1],args=([tk1]),method='lm',jac='2-point',max_nfev=2000)

print(sol)

现在，在文档示例中，我从未见过 objective 函数是按照我的方式编写的。在文档中，objective 函数始终是 returns 残差，而不是残差的平方。此外，在文档中他们从不使用总和！所以我想知道总和和平方是否在 least_squares() 的幕后自动处理？

(2) 也许我的第二个问题是因为我没有理解如何编写 objective 函数。但无论如何，我无法让算法收敛于最小值。我知道这是因为 levenberg 算法是 "greedy" 并在最接近的最小值附近停止，但我认为在不同的初始猜测下我至少能够收敛到大致相同的结果。对最初的猜测稍作改动，我得到了带有不同符号的参数结果。此外，我还没有找到允许算法收敛的初始猜测组合！它总是在找到解决方案之前超时。我什至将函数评估的数量增加到 10,000，看看它是否会。无济于事！

也许有人可以在这里阐明我的错误！我对 python 和 scipy 库还比较陌生！

这是我从此处的视频中自己测量的 tk 的一些样本数据：https://www.youtube.com/watch?v=0Zj_9ypBnzg

tk = [0.52,1.28,2.04,3.17,4.53,6.22]
tk1 = [0.51,1.4,2.09,3,4.42,6.17]
tk2 = [0.63,1.35,2.19,3.02,4.57,6.29]
tk3 = [0.63,1.39,2.23,3.28,4.70,6.32]
tk4 = [0.57,1.4,2.1,3.06,4.53,6.17]

谢谢

Answer 1

1) 是的，正如您所怀疑的，残差的总和和平方是自动处理的。

2) 很难说，因为我对这个问题不是很熟悉（例如，存在多少个局部最小值，什么构成 'reasonable' 结果等）。稍后我可能会调查更多。

但为了好玩，我摆弄了一些值，看看会发生什么。例如，您可以将 1/b 常量替换为独立变量 b_inv，这似乎使结果稳定了很多。这是我用来检查结果的代码。（请注意，为简洁起见，我重写了 objective 函数。它只是利用了 numpy 数组的逐元素操作，而没有改变整体结果。）

import numpy as np
from scipy.optimize import least_squares

def fall_time(k,a,b_inv,c0):
    return (b_inv / a) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))

def parameter_estimation_function(x,tk):
    return np.asarray(tk) - fall_time(k=np.arange(1,len(tk)+1), a=x[0],b_inv=x[1],c0=x[2])

tk_samples = [
    [0.52,1.28,2.04,3.17,4.53,6.22],
    [0.51,1.4,2.09,3,4.42,6.17],
    [0.63,1.35,2.19,3.02,4.57,6.29],
    [0.63,1.39,2.23,3.28,4.70,6.32],
    [0.57,1.4,2.1,3.06,4.53,6.17]
    ]

for i in range(len(tk_samples)):
    sol = least_squares(parameter_estimation_function,[0.1,1.25,-0.1],
        args=(tk_samples[i],),method='lm',jac='2-point',max_nfev=2000)
    print(sol.x)

控制台输出：

[ 0.03621789  0.64201913 -0.12072879]
[  3.59319972e-02   1.17129458e+01  -6.53358716e-03]
[  3.55516005e-02   1.48491493e+01  -5.31098257e-03]
[  3.18068316e-02   1.11828091e+01  -7.75329834e-03]
[  3.43920725e-02   1.25160378e+01  -6.36307506e-03]

SciPy.optimize.least_squares() Objective 函数问题

SciPy.optimize.least_squares() Objective Function Questions

scipy

least-squares

minimize

convergence

levenberg-marquardt