Identify a transfer function from frequency data?

Assume that I want to identify a transfer function from frequency data. Assume that I have a transfer function step response which look like this:

By experience, I would say that this step response is from a second order transfer function on standard form:

$$G(s) = \frac{Y(s)}{U(s)} =\frac{b_0}{a_0s^2 + a_1s + 1 }$$

So I want to find $a_0, a_1, b_0$. Well! I can assume that

$$Y(s) = G(s)U(s)\\ Y(s)(a_1s^2 + a_0s + 1) = b_0U(s) \\ Y(s)a_1s^2 +Y(s) a_0s + Y(s) = b_0U(s)$$

What I know here is the system frequency $s$ and the amplitude input signal $U(s)$ and amplitude output signal $Y(s)$.

I can now use Least Square(LS) to estimate a transfer function by using this:

$$\begin{bmatrix} Y(s_0)\\ Y(s_1)\\ Y(s_2)\\ Y(s_3)\\ \vdots \\ Y(s_n) \end{bmatrix} = \begin{bmatrix} -Y(s_0)s_0^2 & Y(s_0)s_0 & U(s_0) \\ -Y(s_1)s_1^2 & Y(s_1)s_1 & U(s_1) \\ -Y(s_2)s_2^2 & Y(s_2)s_2 & U(s_2) \\ -Y(s_3)s_3^2 & Y(s_3)s_3 & U(s_3) \\ \vdots & \vdots & \vdots \\ -Y(s_n)s_n^2 & Y(s_n)s_n & U(s_n) \\ \end{bmatrix}\begin{bmatrix} a_1 \\ a_0 \\ b_0 \end{bmatrix}$$

Please. Correct me if I'm wrong to this.

Question:

Now I need to find $s$ and $Y(s), U(s)$. One thing is to look at measured data of input $u$ and output $y$.

But how can I do this in MATLAB/Octave? Is there any algorithm or function for this. To find the ampliude and the frequency from measured input and output. The frequency would be the same, but not the amplitude.

In order to identify a second order system of the form

$$G(s) = \dfrac{K\omega_0^2}{s^2+2\zeta\omega_0s+\omega_0^2}$$

You need to determine three values. The first value is the terminal value $f(\infty)=K$ of the time response. Then you determine the peak overshoot $\text{PO}=\dfrac{f_\max(t)-K}{K}$.

By using the formula $$\zeta = \sqrt{\dfrac{\ln^2 \text{PO}}{\pi^2+\ln^2 \text{PO}}}$$ you can determine the parameter $\zeta$.

And the last parameter you need to determine is the time to the first peak $t_\text{P}$. You can use this to determine the natural frequency $\omega_0$ by

$$\omega_0 =\dfrac{\pi}{t_\text{P}\sqrt{1-\zeta^2}}.$$

This procedure should give you quite reasonable results. After having this approximate transfer function you can tweak these parameters.

• How about third order? Or a very large transfer function? It's mutch better to choose the number of poles and zeros, then curve fitting on the transfer function by using least square and data. Dec 25 '17 at 19:16
• How about FFT? Can that help me? Dec 25 '17 at 19:33
• For higher order systems you can determine the explicit solution for the general transfer function by symbolic computation (with Maple, Mathematica, Python or Matlab). Then you simply nonlinear least squares to obtain parameters in the function. Dec 25 '17 at 20:46
• Most of real life system are higher order systems. So I think least square will do an excellent job. But I need to find the amplitude and frequency from sine responses. Dec 25 '17 at 21:05
• I think when it becomes impractical to start approximating the transfer function, you should move to state space techniques. This is just my opinion, not really based on anything. Dec 26 '17 at 15:29

Rewriting the candidate transfer function to the form

$$G(s) = \frac{K}{(s-p_1)(s-p_2)}$$

then the step response of such a system will be

$$y(t) = \frac{K}{p_1\,p_2} + c_1\,e^{p_1\,t} + c_2\,e^{p_2\,t}$$

where $p_1$ and $p_2$ can be complex numbers. So you could perform nonlinear regression in order on the to find ($c_1$, $c_2$), $K$, $p_1$ and $p_2$.

Another option that would work well with any $u(t)$ (assuming sufficient excitation at all frequencies) might be Welch's method. But this will only give you frequency response data and not a transfer function or state space model. However that can still be used to design a controller using loopshaping. But it should be easier to fit a quotient of two polynomials onto the frequency response data then the other time domain data, but this way it should also be easier to estimate the order of the model in general. However if there is noise present in $y(t)$ then you might have to play a little with weights for the data with low signal to noise ratio when fitting.

• Can FFT command in MATLAB be in handy? Dec 25 '17 at 20:43
• @DanielMårtensson Welch's method (pwelch in MATLAB) does use fft as part of its algorithm. Dec 25 '17 at 21:04
• Does it works great to estimate a TF? Dec 26 '17 at 0:16

I'm tired so I say that it cannot be done by estimate a transfer function from frequency response.

The only solution I have found is that if I rewrite this

$$\begin{bmatrix} Y(s_0)\\ Y(s_1)\\ Y(s_2)\\ Y(s_3)\\ \vdots \\ Y(s_n) \end{bmatrix} = \begin{bmatrix} -Y(s_0)s_0^2 & -Y(s_0)s_0 & U(s_0) \\ -Y(s_1)s_1^2 & -Y(s_1)s_1 & U(s_1) \\ -Y(s_2)s_2^2 & -Y(s_2)s_2 & U(s_2) \\ -Y(s_3)s_3^2 & -Y(s_3)s_3 & U(s_3) \\ \vdots & \vdots & \vdots \\ -Y(s_n)s_n^2 & -Y(s_n)s_n & U(s_n) \\ \end{bmatrix}\begin{bmatrix} a_1 \\ a_0 \\ b_0 \end{bmatrix}$$

Into difference equations - They are discrete ODE's:

$$\begin{bmatrix} y(k_0)\\ y(k_1)\\ y(k_2)\\ y(k_3)\\ \vdots \\ y(k_n) \end{bmatrix} = \begin{bmatrix} -y(k_0+2) & -y(k_0+1) & u(k_0) \\ -y(k_1+2) & -y(k_1+1) & u(k_1) \\ -y(k_1+2) & -y(k_2+1) & u(k_2) \\ -y(k_1+2) & -y(k_3+1) & u(k_3) \\ \vdots & \vdots & \vdots \\ -y(k_n+2) & -y(k_n+1) & u(k_n) \\ \end{bmatrix}\begin{bmatrix} a_1 \\ a_0 \\ b_0 \end{bmatrix}$$

Notice that $y(k+2)$ is acceleration, $y(k+1)$ is the derivative and $y(k)$ is the position.

• $y(k)$ should always have the same units, so $y(k+1)$ and $y(k+2)$ can't represent velocity and acceleration, while $y(k)$ represents position. However it can be noted that you can try to fit a discrete transfer function onto data generated from a continues transfer function (but sampled at a fixed rate). A lot of methods related to discrete time model fitting also try to fit some disturbance model as well. Dec 26 '17 at 6:06
• @KwinvanderVeen Thank you. I think that estimate a transfer function is not the best thing for me. I have tried to estimate a state space model from step/impulse/arbitrary responses and it works very good. Except one algorithm, which I going to open a question about now. Right now I using Eigensystem Realization Algorithm and Eigensystem Realization Algorithm Data Correlation and Step-Based Realization. Works really good. Dec 26 '17 at 15:21
• @KwinvanderVeen Hi again! Can you show me a tiny algorithm how to curve fitting. I have implement a code which should do a ARX model. But it works not good at all. Dec 27 '17 at 0:17
• For general nonlinear regression you could just take the sum of the squares of the errors. And then apply gradient descent or quasi-Newton. But if the problem is not convex, them you might only converge to a local minima depending on your initial guess. Dec 27 '17 at 9:10
• @KwinvanderVeen But can I not use this: $$G(s) = \frac{Y(s)}{U(s)}$$ and I know all $s$ and all $G(s)$, all $Y(s)$ and all $U(s)$ which are amplitudes. Then I could just do least square to find $Y(s)G(s) = U(s)$ ? Right? Just choose some polynomial length of $Y(s)$ and $U(s)$? Do I even need to use FFT here if I all ready have measured the sine signal of $u, y$ ? Jan 2 '18 at 22:20