0

我试图通过优化三个未知参数 a、b 和 c0 来最小化高度非线性函数。我正在尝试在 Python 3 中复制赌场轮盘赌球的一些控制方程。

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

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

我的问题是双重的：（1）我正在使用 scipy.optimize.least_squares() 方法='lm'，我不确定如何编写目标函数！现在我写的函数和论文中的完全一样：

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)

现在，在文档示例中，我从未见过以我的方式编写的目标函数。在文档中，目标函数总是返回残差，而不是残差的平方。此外，在文档中，他们从不使用总和！所以我想知道总和和平方是否在引擎盖下自动处理least_squares()？

(2) 也许我的第二个问题是由于我没有理解如何编写目标函数。但无论如何，我无法让算法收敛到最小值。我知道这是因为 levenberg 算法是“贪婪的”并且在最接近的最小值附近停止，但我认为考虑到不同的初始猜测，我至少能够收敛到大致相同的结果。随着最初猜测的轻微变化，我得到了不同符号的参数结果。此外，我还没有找到允许算法收敛的初始猜测组合！它总是在找到解决方案之前超时。我什至将函数评估的数量增加到 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]

谢谢

4

1 回答 1

2

1）是的，因为您怀疑残差的总和和平方是自动处理的。

2）很难说，因为我对这个问题并不十分熟悉（例如，存在多少个局部最小值，什么构成了“合理”的结果等）。稍后我可能会进行更多调查。

但是为了好玩，我摆弄了一些值来看看会发生什么。例如，您可以将1/b常量替换为独立变量b_inv，这似乎使结果稳定了很多。这是我用来检查结果的代码。（请注意，为简洁起见，我重写了目标函数。它只是利用了 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]
于 2017-12-26T14:32:52.833 回答