Fitting transfer function models in Scipy.signal

Question

I am using curve_fit to fit a step response of a first order dynamic system to estimate the gain and time constant. I use two approaches. First approach is to fit the curve generated from the function , in the time domain .

# define the first order dynamics in  the time domain
def model(t,gain,tau):
    return (gain*(1-exp(-t/tau)))

#define the time intervals
time_interval = linspace(1,100,100)

#genearte the output using the model with gain= 10 and tau= 4 
output= model(t,10,4)

# fit to output and estimate  parameters - gain and tau
par = curve_fit(time_interval, output)

Now checking par reveals an array of 10 and 4 which is perfect.

The second approach is to estimate gain and time constant by fitting to a step response of a LTI system The LTI System is defined as a transfer function with numerator and denominator.

#define function as a step response of a LTI system .
# The argument x has no significance here,
# I have included because , the curve_fit requires passing "x" data to the function

def model1(x ,gain1,tau1):
    return lti(gain1,[tau1,1]).step()[1]

#generate output using the above model
output1 = model1(0,10,4)

par1 = curve_fit(model1,1,output1)

now checking par1 reveals an array of [ 1.00024827, 0.01071004] which is wrong. What is wrong with my second approach here? Is there more efficient way of estimating the transfer function coefficients from the data by curve_fit

Thank you

Answer 1

The first three arguments to curve_fit are the function to be fit, the xdata and the ydata. You have passed xdata=1. Instead you should give it the time values associated with output1.

One way to do that is to actually use the first argument in the function model1, like you did in model(). For example:

import numpy as np
from scipy.signal import lti
from scipy.optimize import curve_fit


def model1(x, gain1, tau1):
    y = lti(gain1, [tau1, 1]).step(T=x)[1]
    return y

time_interval = np.linspace(1,100,100)

output1 = model1(time_interval, 10, 4)

par1 = curve_fit(model1, time_interval, output1)

I get [10., 4.] for the parameters, as expected.