Search code examples
pythonnumpyscipycurve-fitting

Curve Fitting with SciPy failing to give a correct fit

I have two NumPy arrays x_data and y_data. When I try to fit my data using a second order step response model and scipy.optimize.curve_fit with this code:

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}")

This is the equation I am using to fit the data enter image description here

I get the following fit. The K_fitted(gain) and zeta_fitted(dampening coefficient) are within realistic values but omega_n_fitted is way to large. enter image description here

What is the problem?

Edit1: Adding image over underdamped and overdamped systems for context. I am fitting the underdamped system enter image description here

Solution

  • COMMENT :

    The model equation appears not well compatible with the data. The fitting is far to be correct.

    A different model equation is considered below. The variable t is repaced by a variable x=ln(t-t0). Of course this isn't significant on physical viewpoint. This is purely mathematical.

    enter image description here