Question

我有两个 NumPy 数组x_data和y_data.当我尝试使用二阶阶跃响应模型和scipy.optimize.curve_fit以下代码来拟合数据时：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

x_data = np.array([51.5056, 51.5058, 51.5061, 51.5064, 51.5067, 51.5069, 51.5072,
       51.5075, 51.5078, 51.5081, 51.5083, 51.5086, 51.5089, 51.5092,
       51.5094, 51.5097, 51.51  , 51.5103, 51.5106, 51.5108, 51.5111,
       51.5114, 51.5117, 51.5119, 51.5122, 51.5125, 51.5128, 51.5131,
       51.5133, 51.5136, 51.5139, 51.5142, 51.5144, 51.5147, 51.515 ,
       51.5153, 51.5156, 51.5158, 51.5161, 51.5164, 51.5167, 51.5169,
       51.5172, 51.5175, 51.5178, 51.5181, 51.5183, 51.5186, 51.5189,
       51.5192, 51.5194, 51.5197, 51.52  , 51.5203, 51.5206, 51.5208,
       51.5211, 51.5214, 51.5217, 51.5219, 51.5222, 51.5225, 51.5228,
       51.5231, 51.5233, 51.5236, 51.5239, 51.5242, 51.5244, 51.5247,
       51.525 , 51.5253, 51.5256, 51.5258, 51.5261, 51.5264, 51.5267,
       51.5269, 51.5272, 51.5275, 51.5278, 51.5281, 51.5283, 51.5286,
       51.5289, 51.5292, 51.5294, 51.5297, 51.53  , 51.5303, 51.5306,
       51.5308, 51.5311, 51.5314, 51.5317, 51.5319, 51.5322, 51.5325,
       51.5328, 51.5331, 51.5333, 51.5336, 51.5339, 51.5342])

y_data = np.array([2.99 , 2.998, 3.024, 3.036, 3.038, 3.034, 3.03 , 3.025, 3.02 ,
       3.016, 3.012, 3.006, 3.003, 3.   , 2.997, 2.995, 2.993, 2.99 ,
       2.989, 2.987, 2.986, 2.985, 2.983, 2.983, 2.982, 2.98 , 2.98 ,
       2.979, 2.978, 2.978, 2.976, 2.977, 2.976, 2.975, 2.975, 2.975,
       2.975, 2.974, 2.975, 2.974, 2.974, 2.974, 2.973, 2.974, 2.973,
       2.974, 2.973, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974,
       2.973, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974,
       2.975, 2.974, 2.975, 2.975, 2.976, 2.976, 2.976, 2.976, 2.975,
       2.976, 2.976, 2.977, 2.977, 2.976, 2.977, 2.977, 2.977, 2.977,
       2.978, 2.978, 2.978, 2.979, 2.978, 2.978, 2.979, 2.979, 2.979,
       2.979, 2.979, 2.98 , 2.98 , 2.98 , 2.98 , 2.98 , 2.981, 2.981,
       2.981, 2.981, 2.982, 2.981, 2.982])

# Second order function definition
def second_order_step_response(t, K, zeta, omega_n):
    omega_d = omega_n * np.sqrt(1 - zeta**2)
    phi = np.arccos(zeta)
    return K * (1 - (1 / np.sqrt(1 - zeta**2)) * np.exp(-zeta * omega_n * t) * np.sin(omega_d * t + phi))

# Initial Guess of parameters 
K_guess = y_data.max() - y_data.min()
zeta_guess = 0.5  # Typically between 0 and 1 for underdamped systems
omega_n_guess = 2 * np.pi / (x_data[1] - x_data[0])  # A rough estimate based on data sampling rate

# Fit the model with increased maxfev and parameter bounds
params, covariance = curve_fit(
    second_order_step_response, x_data, y_data,
    p0=[K_guess, zeta_guess, omega_n_guess],
    maxfev=5000,  # Increase max function evaluations
    bounds=([0, 0, 0], [np.inf, 1, np.inf])  # Bounds for K, zeta, and omega_n
)
K_fitted, zeta_fitted, omega_n_fitted = params

# Generate data for the fitted model
x_fit = np.linspace(min(x_data), max(x_data), 300)  # More points for a smoother line
y_fit = second_order_step_response(x_fit, K_fitted, zeta_fitted, omega_n_fitted)

# Plot
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, color='red', label='Original Data')  # Original data
plt.plot(x_fit, y_fit, 'blue', label='Fitted Model')  # Fitted model
plt.title('Fitting Second Order System Step Response to Data')
plt.xlabel('Time')
plt.ylabel('Y')
plt.legend()
plt.show()

# Display fitted parameters
print(f"Fitted Parameters: Gain (K) = {K_fitted}, Damping Ratio (zeta) = {zeta_fitted}, Natural Frequency (omega_n) = {omega_n_fitted}")

这是我用来拟合数据的方程

我得到以下适合。（增益K_fitted）和zeta_fitted（阻尼系数）在实际值范围内，但omega_n_fitted太大了。

问题是什么？

Edit1：在欠阻尼和过阻尼系统上添加图像以获取上下文。我正在安装欠阻尼系统

您对参数的最初猜测与正确的解决方案相差太远，优化就会变得混乱。它应该警告您适配不良/失败。时间轴上的加法常数 51.5 也无助于数值稳定性。从峰的半高宽推测参数，它可能有机会。 — 
你的方程是从哪里来的？为什么定义了 aphi但不使用它？ — 
@jared我正在重写方程并忘记删除phi.我编辑了代码以包含phi.我添加了一张描述方程的图片 — 
我可以得到一个合适的（使用修改后的方程），但由于那条弯曲的尾巴，它不是那么好。 — 
您的数据看起来频率随时间变化，而您的模型具有恒定频率。 —

Accepted Answer

评论：

模型方程似乎与数据不太兼容。拟合远不正确。

下面考虑一个不同的模型方程。变量 t 由变量 x=ln(t-t0) 替换。当然，从物理角度来看，这并不重要。这纯粹是数学问题。

3

这是对数据的巧妙的实证拟合。我的拟合首选 51.5056，并且使用经典的二阶阻尼模型无法正确拟合尾部。

–
@JJAcquelin。哇，真是个巫师。我使用的模型确实似乎不适合数据。您的模型能够捕获数据的尾部和剩余部分。我将尝试从物理角度理解这些方程为何有效以及它们在该系统中所描述的内容。谢谢你！

–
我怀疑经验模型是否具有物理意义。也许人们可以找到其他也适合数据的经验模型，但无法说出哪一个从物理角度来看更好。我会更有信心通过研究所涉及的现象来建立一个模型。祝你好运！

–

|

Answer 2

下面给出了一个相当标准的物理模型。

对于现有数据（但请参阅本文底部），最佳拟合似乎是临界阻尼（阻尼比 zeta=1）。尝试与欠阻尼或过阻尼解决方案进行最佳配合会得到 zeta–>1，因此您不妨直接采用临界阻尼。

线性阻尼模型（机械或电气）的平衡位移解供参考。自由常数 A、B 来自边界条件（即，去除力的点）。

过阻尼和临界阻尼的平衡 y 值为 3.01，omega=76.47。这些解决方案在图表上看起来几乎无法区分。

欠阻尼解决方案可以获得初始“凹凸”，但无法跟踪尾部。临界阻尼和过阻尼解决方案获得大量数据，但不是开始（可能仍包括一些初始强迫）。

我从 x 数据中减去了初始时间 t0。没有真正需要适应这一点：如果将其作为自由参数保留，那么它只会被吸收到乘以指数项的常数中或改变正弦项中的相位。

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

x_data = np.array([51.5056, 51.5058, 51.5061, 51.5064, 51.5067, 51.5069, 51.5072,
       51.5075, 51.5078, 51.5081, 51.5083, 51.5086, 51.5089, 51.5092,
       51.5094, 51.5097, 51.51  , 51.5103, 51.5106, 51.5108, 51.5111,
       51.5114, 51.5117, 51.5119, 51.5122, 51.5125, 51.5128, 51.5131,
       51.5133, 51.5136, 51.5139, 51.5142, 51.5144, 51.5147, 51.515 ,
       51.5153, 51.5156, 51.5158, 51.5161, 51.5164, 51.5167, 51.5169,
       51.5172, 51.5175, 51.5178, 51.5181, 51.5183, 51.5186, 51.5189,
       51.5192, 51.5194, 51.5197, 51.52  , 51.5203, 51.5206, 51.5208,
       51.5211, 51.5214, 51.5217, 51.5219, 51.5222, 51.5225, 51.5228,
       51.5231, 51.5233, 51.5236, 51.5239, 51.5242, 51.5244, 51.5247,
       51.525 , 51.5253, 51.5256, 51.5258, 51.5261, 51.5264, 51.5267,
       51.5269, 51.5272, 51.5275, 51.5278, 51.5281, 51.5283, 51.5286,
       51.5289, 51.5292, 51.5294, 51.5297, 51.53  , 51.5303, 51.5306,
       51.5308, 51.5311, 51.5314, 51.5317, 51.5319, 51.5322, 51.5325,
       51.5328, 51.5331, 51.5333, 51.5336, 51.5339, 51.5342])

y_data = np.array([2.99 , 2.998, 3.024, 3.036, 3.038, 3.034, 3.03 , 3.025, 3.02 ,
       3.016, 3.012, 3.006, 3.003, 3.   , 2.997, 2.995, 2.993, 2.99 ,
       2.989, 2.987, 2.986, 2.985, 2.983, 2.983, 2.982, 2.98 , 2.98 ,
       2.979, 2.978, 2.978, 2.976, 2.977, 2.976, 2.975, 2.975, 2.975,
       2.975, 2.974, 2.975, 2.974, 2.974, 2.974, 2.973, 2.974, 2.973,
       2.974, 2.973, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974,
       2.973, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974, 2.974,
       2.975, 2.974, 2.975, 2.975, 2.976, 2.976, 2.976, 2.976, 2.975,
       2.976, 2.976, 2.977, 2.977, 2.976, 2.977, 2.977, 2.977, 2.977,
       2.978, 2.978, 2.978, 2.979, 2.978, 2.978, 2.979, 2.979, 2.979,
       2.979, 2.979, 2.98 , 2.98 , 2.98 , 2.98 , 2.98 , 2.981, 2.981,
       2.981, 2.981, 2.982, 2.981, 2.982])

x_data = x_data - x_data[0]


# Second order function definition
def underDamped(t, y0, A, B, omega, zeta ):
    omega_d = omega * np.sqrt( 1 - zeta ** 2 )
    return y0 + A * np.exp( -zeta * omega * t ) * np.sin( omega_d * t + B )
under0 = [ 3, 0.1, 0.0, 50.0, 0.1 ]
punder, cv = curve_fit( underDamped, x_data, y_data, p0=under0, bounds=([0,-np.inf,-np.inf,0.0,0.0],[10,np.inf,np.inf,np.inf,1.0]), maxfev=10000 )

def overDamped(t, y0, A, B, omega, zeta ):
    lambda1 = omega * ( zeta - np.sqrt( zeta ** 2 - 1 ) )
    lambda2 = omega * ( zeta + np.sqrt( zeta ** 2 - 1 ) )
    return y0 + A * np.exp( -lambda1 * t ) + B * np.exp( -lambda2 * t )
over0 = [ 3, 0.1, 0.1, 50.0, 2.0 ]
pover, cv = curve_fit( overDamped, x_data, y_data, p0=over0, maxfev=10000 )

def criticallyDamped(t, y0, A, B, omega ):
    return y0 + ( A + B * t ) * np.exp( -omega * t )
crit0 = [ 3, 0.1, 0.1, 50.0 ]
pcrit, cv = curve_fit( criticallyDamped, x_data, y_data, p0=crit0, maxfev=10000 )

# Generate data for the fitted model
x_fit = np.linspace(min(x_data), max(x_data), 300)  # More points for a smoother line
u_fit = underDamped     ( x_fit, *punder )
o_fit = overDamped      ( x_fit, *pover  )
c_fit = criticallyDamped( x_fit, *pcrit  )

# Plot
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, color='black', label='Original Data')
plt.plot(x_fit, u_fit, 'red',   label='Underdamped')
plt.plot(x_fit, o_fit, 'blue' , label='Overdamped')
plt.plot(x_fit, c_fit, 'green', label='Critically damped')
plt.title('Fitting Second Order System Step Response to Data')
plt.xlabel('Time')
plt.ylabel('Y')
plt.legend()
plt.show()

# Display fitted parameters
print( "Underdamped:       y0, A, B, omega, zeta = ", *punder )
print( "Overdamped:        y0, A, B, omega, zeta = ", *pover  )
print( "Critically-damped: y0, A, B, omega = "      , *pcrit  )

输出：

Underdamped:       y0, A, B, omega, zeta =  2.976449192328566 -45.10851382377807 -3.1413665103191675 815.1575670585344 0.9999946970217488
Overdamped:        y0, A, B, omega, zeta =  3.011323385679147 -4.798576387293144 4.8163356019586825 76.46587506775039 1.000079837533067
Critically-damped: y0, A, B, omega =  3.011320122667347 0.017762157049911867 -9.290153164255033 76.47022695793991

编辑：根据 Martin Brown 的评论，我也尝试删除前 4 点。临界阻尼（阻尼比 zeta=1）非常适合这种情况。平衡位置 (y0=2.982) 和自然圆频率 (omega=172.686) 的修正值。

Critically-damped: y0, A, B, omega =  2.982416646620049 0.055681833492928545 -10.697532138927013 172.6861383346045

我怀疑这些数据的前4点是在施加一些外力的过程中的，因此不属于解决方案。假装它是临界阻尼的，可以让拟合的上升时间稍微陡一些。我得到了类似的拟合，无论我做什么，我都无法在出口翼上重现线性上升的基线。 OP 没有告诉我们有关他的实验的一些细节。拟合 L2 与 L1 范数也没有帮助（后者将残差大部分移动到基线上）。 — 
“我怀疑这些数据的前 4 点是在施加某些外力时的”——是的，我同意这一点。我仍然认为随后的行为（在没有明确的振荡频率的情况下）意味着临界阻尼或过阻尼。如果删除前 4 个点，则临界阻尼非常适合该行为。 — 
@马丁布朗。该数据是瞬态电压响应的摘录。简而言之，我们从低电流到高电流并测量电压与时间的关系。你是对的。在前 4 个数据点之前约 1 秒施加外力（高电流）。这里显示的所有数据都是在较高电流下的。我们得到的电压与时间数据是系统响应 — 
@impedance_gatto，电路中有什么？如果您知道电阻、电容和电感，您就可以计算出固有频率和阻尼比。就目前情况而言，看起来有人已经完美地将其设置为临界阻尼（这通常是一个很好的结果）。 — 
后来的配合看起来更有说服力。临界阻尼是运气好还是设计使然？您始终可以为其加载一些额外的电容或电感作为最终测试。 —

python – 使用SciPy进行曲线拟合未能给出正确的拟合 – VoidCC

2 个回答
2

2 个回答 2

2 个回答
2