4. Fitting step responses

It is often prohibitively expensive to develop first principle models of processes and therefore it is very common to estimate low order transfer functions directly from plant data. This is simple to do if we have access to step test results.

import control
import numpy
import matplotlib.pyplot as plt
import tbcontrol
tbcontrol.expectversion("0.1.10")
%matplotlib inline

Let’s start with a higher order process to generate our “real data”

Greal = control.tf([1, 2], [2, 3, 4, 1])

ts, ys = control.step_response(Greal)

Remember that the real data will not necessarily start at zero, so we’ll add in some value for the initial output. We will also add some normally distributed noise to represent measurement error.

yinitial = 10
measurement_noise = numpy.random.randn(len(ys))*0.05

ym = ys + yinitial + measurement_noise

plt.scatter(ts, ym)
plt.plot(ts, ys + yinitial, color='red')

[<matplotlib.lines.Line2D at 0x1c1f8f10b8>]

import scipy.optimize

We’ll fit a first order plus dead time and second order plus dead time model. The tbcontrol.responses library contains the analytic formulae for these step responses.

from tbcontrol.responses import fopdt, sopdt

It is very important to have a good idea of the initial parameter values. Interaction makes it easy to figure out.

from ipywidgets import interact

def resultplot(K, tau, theta, y0):
    plt.scatter(ts, ym)
    plt.plot(ts, fopdt(ts, K, tau, theta, y0), color='red')
    plt.show()

interact(resultplot,
         K=(1., 10.),
         tau=(0., 10.),
         theta=(0., 10.),
         y0=(0., 20.));

We can use the scipy.optimize.curve_fit tool to do this fit just like when we did regression without time.

[K, tau, theta, y0], _ = scipy.optimize.curve_fit(fopdt, ts, ym, [2, 4, 1, 10])
[K, tau, theta, y0]

[1.954094509563327, 2.8181973943514884, 0.6113717655974688, 10.031145239708668]

The parameters for the second order model should be similar, with a smaller time constant and overdamped nature.

[K_2, tau_2, zeta_2, theta_2, y0_2], _ = scipy.optimize.curve_fit(sopdt, ts, ym, [2, 2, 1.5, 1, 10])
[K_2, tau_2, zeta_2, theta_2, y0_2]

[1.9798056401724866,
 0.8305183954384022,
 1.8452576915795178,
 0.3199080792764559,
 10.005597510469169]

plt.figure(figsize=(10, 5))
plt.scatter(ts, ym, label='Data')
plt.plot(ts, fopdt(ts, K, tau, theta, y0), color='red', label='FOPDT fit')
plt.plot(ts, sopdt(ts, K_2, tau_2, zeta_2, theta_2, y0_2), color='red', label='SOPDT fit')
plt.plot(ts, ys + 10, color='blue', label='Original')
plt.legend(loc='best')

<matplotlib.legend.Legend at 0x1c1fd1c710>