Search code examples
pythonscipycurve-fittinglti

Fit First Order Transfer Function (Curve Fitting)

I have the following code to fit the curve and show me the time constant tau. However, when I try to fit the experimental data it doesn't work. Anybody can help?

from scipy.signal import lti
import pylab as plt
from scipy.optimize import curve_fit
import math as m

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

def get_time_constants(time_vector, power_vector):
    time_constant = []
    for i in range(len(power_vector)):
        # fit to output and estimate  parameters - gain and tau
        par1 = curve_fit(model, time_vector[i], power_vector[i])
        y_fit = model(time_vector[i], par1[0][0], par1[0][1])

        time_constant.append(par1[0][1])
        time_constant_mean = sum(time_constant) / len(time_constant)
        plt.plot(time_vector[i], power_vector[i])
        plt.plot(time_vector[i], y_fit, label='Time Constant: %.2f s' %par1[0][1])
        plt.show()

power = [[0.0, 68.13, 108.445, 133.43, 158.56, 164.575, 168.26, 172.035, 172.94, 173.795, 173.96, 174.145, 174.165, 174.195, 174.215, 174.29, 174.305, 174.325]]
time = [[0.0, 1.0, 2.0, 3.0999999999999943, 5.099999999999994, 6.099999999999994, 7.099999999999994, 9.099999999999994, 10.199999999999989, 12.199999999999989, 13.199999999999989, 14.199999999999989, 16.19999999999999, 17.299999999999983, 19.299999999999983, 20.30000000000001, 21.30000000000001, 23.30000000000001]]

get_time_constants(time, power)

This is the result I am getting: enter image description here

Solution

  • As so many times ... just a question of starting values

    def get_time_constants(time_vector, power_vector):
    time_constant = []
    for t, p  in zip( time_vector, power_vector ):
        guess_amp = max( p )
        guess_tau = guess_amp / ( ( p[1] - p[0]) / ( t[1] - t[0] ) )
        par1, _ = curve_fit(model, t, p, p0=[ guess_amp, guess_tau ])
        y_fit = model(t, *par1)
        time_constant.append(par1[1])
        time_constant_mean = sum(time_constant) / len(time_constant)
        plt.plot( t, p, ls='', marker='+' )
        plt.plot( t, y_fit, label='Time Constant: %.2f s' %par1[1])
        plt.legend(loc=0)
        plt.show()

    works just fine.