scipy-curve-fit

Using scipy.optimize.curve_fit for Parameter Estimation

Overview

scipy.optimize.curve_fit is a tool for fitting models to experimental data using nonlinear least squares optimization.

Basic Usage

from scipy.optimize import curve_fit
import numpy as np

# Define your model function
def model(x, param1, param2):
    return param1 * (1 - np.exp(-x / param2))

# Fit to data
popt, pcov = curve_fit(model, x_data, y_data)

# popt contains the optimal parameters [param1, param2]
# pcov contains the covariance matrix

Fitting a First-Order Step Response

import numpy as np
from scipy.optimize import curve_fit

# Known values from experiment
y_initial = ...  # Initial output value
u = ...          # Input magnitude during step test

# Define the step response model
def step_response(t, K, tau):
    """First-order step response with fixed initial value and input."""
    return y_initial + K * u * (1 - np.exp(-t / tau))

# Your experimental data
t_data = np.array([...])  # Time points
y_data = np.array([...])  # Output readings

# Perform the fit
popt, pcov = curve_fit(
    step_response,
    t_data,
    y_data,
    p0=[K_guess, tau_guess],      # Initial guesses
    bounds=([K_min, tau_min], [K_max, tau_max])  # Parameter bounds
)

K_estimated, tau_estimated = popt

Setting Initial Guesses (p0)

Good initial guesses speed up convergence:

# Estimate K from steady-state data
K_guess = (y_data[-1] - y_initial) / u

# Estimate tau from 63.2% rise time
y_63 = y_initial + 0.632 * (y_data[-1] - y_initial)
idx_63 = np.argmin(np.abs(y_data - y_63))
tau_guess = t_data[idx_63]

p0 = [K_guess, tau_guess]

Setting Parameter Bounds

Bounds prevent physically impossible solutions:

bounds = (
    [lower_K, lower_tau],    # Lower bounds
    [upper_K, upper_tau]     # Upper bounds
)

Calculating Fit Quality

R-squared (Coefficient of Determination)

# Predicted values from fitted model
y_predicted = step_response(t_data, K_estimated, tau_estimated)

# Calculate R-squared
ss_residuals = np.sum((y_data - y_predicted) ** 2)
ss_total = np.sum((y_data - np.mean(y_data)) ** 2)
r_squared = 1 - (ss_residuals / ss_total)