SciPy.optimize.least_squares()目标函数问题

问题描述 投票:0回答:1

我试图通过优化三个未知参数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]

谢谢

scipy least-squares minimize convergence levenberg-marquardt
1个回答
1
投票

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]
© www.soinside.com 2019 - 2024. All rights reserved.