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

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：在欠阻尼和过阻尼系统上添加图像以获取上下文。我正在安装欠阻尼系统

Answer 1

实际上，最佳拟合似乎是临界阻尼（zeta=1）。即使尝试使用过阻尼解决方案也会得到 zeta-->1。

过阻尼和临界阻尼的平衡 y 值为 3.01，omega=76.47。

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 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, 1.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, 1.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
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='red', label='Original Data')
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( "Over-damped: y0, A, B, omega, zeta = ", *pover )
print( "Critically-damped: y0, A, B, omega = ", *pcrit )

输出：

Over-damped: y0, A, B, omega, zeta =  3.0113230036207916 -4.979959895551821 4.997719482218953 76.46631309103108 1.0000741371536719
Critically-damped: y0, A, B, omega =  3.0113157317741974 0.017766928030121528 -9.289980190715237 76.47399440858744