# How can i extract transfer function an unknown nonlinear system?

I am trying to do PID control for my electroservo motor system by using nichols ziegler tuning method. My system has SSI encoder output for motor feedback mechanism. I will use this knowledge for control. According to nichols ziegler method i must know transfer function of my system. But i can not find its equation exactly. So how can i extract its transfer function? I need a methodology for this. Can i extract T.F. by using Matlab/Simulink or LAbview?

• No - all the simulation tools need the transfer function as an INPUT. But there are many technical cobtributions and papers dealing with dc servo motors and the corresponding tranfer functions (google: transfer function for dc motor). – LvW Nov 26 '15 at 9:22
• Sir, while searching i saw "tfest" function on matlab. But for using this function, i must know poles and zeros of my system. I can look similar system to identify poles and zeros. I guess there is a method for extracting T.F. First step, i can supply step function to my system input and i observe and record output data. Both input and output can transform to "s domain" by using laplace transform. And tf(s)=Output(s)/Input(s). Is that way correct? – Cem Nov 26 '15 at 11:46
• The Ziegler-Nichols method is the only method where you don't have to know the transfer function for setting the parameters, that is why this method exists. – Marko Buršič Nov 26 '15 at 14:46
• You should tune and then run chirp test to extract controller performance. That will tell you how well your controller is able to follow the setpoint across the whole frequency range given the specific controller gains. That is a hard metric that is very easy to use to compare different gains and to pick the best tuning that keeps input to output gain at zero for as high frequency as possible (ie perfect setpoint following). – Martin Feb 4 '18 at 14:22
• You can use functions I have provided below (freq_response_frd) to generate a frequency response from a chirp input and measured output and then plot it using bode(freq_response_frd(...)) to get the bode plot of the system response. – Martin Feb 4 '18 at 14:22

The concept of Transfer Function is only defined for linear time invariant systems. Nonlinear system models rather stick to time domain descriptions as nonlinear differential equations rather than frequency domain descriptions.

But in terms of current-in, speed out, your motor-encoder system is close enough to a linear system that you really don't need to concern yourself with nonlinear aspects (unless you are trying to control shaft angle to micro-radian precision!).

Perhaps the easiest way to obtain a linear model is to apply a simple proportional feedback control tuned to get the loop stable, then record input-output data to a step response. Then fit the data to the closed loop transfer function. From the closed loop transfer function you can calculate the open loop transfer function, factor out the proportional gain and voila - your motor model! A simple linear DC motor model looks like: $$\frac{\omega}{i}=\frac{K_T}{Js+B}$$ where $$K_T$$ is the torque constant of the motor, $$J$$ is the motro shaft and load inertia and $$B$$ is the linear viscous damping of the motor bearings

Perhaps your motor supplier already specifies these parameters in which case you don't have to test - you can write the model directly.

Note that even if you are using a permanent magnet synchronous motor, in feedback with a stiff current controller, the model approaches the model of the DC (brush) motor.

• Thank you for your reply, but what do you mean by saying " Then fit the data to the closed loop transfer function"? What is mathematical response of this sentence? – Cem Dec 4 '15 at 14:37
• @Cem 'Fitting' simply means determining by some method what choice of KT, J and B provide some minimal difference between the model's response to input data and the actual system's output response. Methods could include such things as least squares, nonlinear optimization, gradient search, etc. – docscience Dec 4 '15 at 20:54

Ziegler-Nicholls tuning does not require the TF to be known - that's the whole point of the method.

• Yes, Z-N is an interactive tuning method that works if the system is similar enough to the assumed model. No prior knowledge of system parameters is required. – Spehro Pefhany Nov 26 '15 at 16:00

You can definitely try System Identification Toolbox of Matlab. Official page says

You can use time-domain and frequency-domain input-output data to identify continuous-time and discrete-time transfer functions, process models, and state-space models.

Which is what you are looking for.

• Can you send me a message? We can better communicate with our native language. I am from Ankara ;-) – Cem Nov 26 '15 at 21:00
• Look at my reply. I have attached octave code you can use. – Martin Feb 4 '18 at 14:17

The first, like @LvW suggested, is to use the motor specifications such as in the illustration. This may be bore accurate but also more difficult because not all the specks will be available.

The second involves what you suggested. By applying a voltage to the motor and recording how its speed behaves. In effect this creates a step response that can be analyzed to find the specs for either a first order or second order system. If the system oscillates then a second order is what you want. If the system does not oscillate then you may be able to use a first order.

• That was perfect explanation! Thank you for this. I am also curious about how do i know its order? What if it is fourth order? By the way the system that i mentioned is "shaking table" project. – Cem Nov 26 '15 at 20:26
• @Cem The first solution will give you the correct order system. With the second solution there are limits. If the system is overdamped then even a third or fourth order system will look like a first order with the dominant pole driving the response. Normally we try to fit a system to a second order just because it's relatively easy to figure out the natural frequency and the dampening factor. Just be sure to have enough phase and gain margin so that if the model is not 100% accurate the controller can tolerate it. – vini_i Nov 26 '15 at 20:46

There are two very good methods for estimating transfer functions. Look up moen4 and fitfrd.

To use moen4 you need basically input and an output of a test. The algorithm then computes the transfer function that best fits the data. The results tend to be pretty good for some systems, less so for systems that have significant non-linear behavior (for which a linear transfer function does not exist).

Here is code that you can use for both frd fit and moen4 fit. You can plot the data of freq_response_frd (frd object) directly using the bode() function to get a bode plot of your input data. Your input data must have sufficient frequency coverage so use a chirp signal that increases with frequency in time and collect the resulting response in another array. Then pass both arrays into the id_model_moen and you will get your transfer function back.

I typically limit the frequencies that are analyzed because if you plot the full range returned by fft you will get a lot of noise outside of the range for which you even have test data - so that's useless part of results.

function [mag, phase, f] = freq_response_mag_phase(out, in, t, freqlim)
dt = (t(end) - t(1)) / (length(t) - 1);
NFFT = length(t);
Fs = 1.0 / dt;
fb = fft(out, NFFT);
fa = fft(in, NFFT);
f = [0:NFFT-1]*Fs/NFFT * 2 * pi;
% find first bin after our test range. We will discard bins after it.
ix = ceil(NFFT/2);
if(exist("freqlim", "var"))
ix = find(f>freqlim,1);
end
f = f(1:ix);
mag = abs(fb(1:ix)) ./ abs(fa(1:ix));
phase = unwrap(angle(fb(1:ix))) - unwrap(angle(fa(1:ix)));
end

function response = freq_response_frd(out, in, t, freqlim)
if(exist("freqlim", "var"))
[mag, phase, f] = freq_response_mag_phase(out, in, t, freqlim);
else
[mag, phase, f] = freq_response_mag_phase(out, in, t);
end
response = frd(mag .* exp(1i .* phase), f);
end

function sys_tf = id_model_frd(out, in, t, nr)
resp = freq_response_frd(out, in, t);
sys = fitfrd(resp, nr);
[b, a] = ss2tf(minreal(sys));
sys_tf = tf(b, a);
end

function sys_tf = id_model_moen(out, in, t, nr)
dt = (t(end) - t(1)) / (length(t) - 1);
sys = moen4(iddata(out, in, dt), nr);
[b,a] = ss2tf(d2c(sys));
sys_tf = tf(b, a);
end