Estimation of transfer function

Hello,

I would like to estimate a transfer function based on the experimental data (magnitude, phase and frequency. I was trying to use pem function but it does not work well. For example in my system I have a derivative and my transfer function should be something like: (as^2+bs)/(a1s^3+b1s^2+c1s+d1). I cannot get it in the estimation model. It gives me something like as^2+bs+c. I do not want "c" parameter. There is any estimation function that allows me to settle my searching function not only number of orders? I mean user0searching parameter? If is there any how I can do it?

Thank you for any response.

Kamil

Answer 1

[Arnaud Miege]

For frequency-domain data, there are three types of id models supported by the System Identification Toolbox(see documentation):

• Low-order transfer functions
• Input-output polynomial models
• State-space models

For your application, I would recommend the latter. Choose the correct order for your system and once you have your estimated state-space matrices, you can use the tf command from the Control System Toolbox to convert it to a transfer function.

HTH,

Arnaud

Answer 2

[Rajiv Singh]

There is also the "grey box" approach that lets you parametrize your model any way you like. For this, write a MATLAB function that takes in a, b, a1, b1, c1 and d1 as input arguments and returns a state-space form of the model as output arguments. Use this file to create an "idgrey" model object whose parameters (your a, b, a1,..) can be estimated to fit data. You will need to derive the expressions for SS matrices to begin with. A controls textbook (e.g., Linear Analysis by Kailath) would show you some common forms.

See: http://www.mathworks.com/help/toolbox/ident/ug/bqs6lf8.html

for details on grey-box modeling. Also see http://www.mathworks.com/products/sysid/demos.html?file=/products/demos/shipping/ident/iddemo7.html

Answer 3

[Kamil]

Hello,

Thank you both of you for answers. I had much work and I could not answer at all. I will try to do it now. The links are to Matlab 2011 I have Matlab 2007b. I hope it will work as well. I will give you answer if I get what I want.

Kamil

Answer 4

[Kamil]

Hello again,

I did test of the estimation but I cannot get the right answer. Well, I can but it strictly depends on the initial points. As I know the function I know what to put. Could you check if I am going in the right direction?

G=tf([1 1],[1 1 1]);

[num,den] = tfdata(G,'v');

[a,b,c,d]=tf2ss(num,den);

[m,ph,f]=bode(G);
zfr = m(:).*exp(i*ph(:)*pi/180);
Ts = 0;
data = idfrd(zfr,f,Ts);

A = [1 1; 1 1];
B = [1; 1];
C = [1 1];
D = [0];

m = idss(A,B,C,D,'Ts',0);

m.As = [NaN NaN; NaN NaN];
m.Bs = [NaN; NaN];
m.Cs = [NaN NaN];
m.Ds = [NaN];

res = pem(data,m);

A=res.a;
B=res.b;
C=res.c;
D=res.d;

[num1, den1] = ss2tf(A, B, C, D);

tf(num1,den1)

Thank you.

Kamil

Answer 5

[Kamil]

Hello again,

I changed a file to my own measured data but unfortunately Matlab has problem to find something similar. Anybody has any idea how I can fix it?

function estimation_data

format short e

f=[6.2832e-001  1.2566e+000  1.8850e+000  2.5133e+000  3.1416e+000  3.7699e+000  4.3982e+000  5.0265e+000  5.6549e+000  6.2832e+000  6.9115e+000...
    7.5398e+000  8.1681e+000  8.7965e+000  9.4248e+000  1.0053e+001  1.0681e+001  1.1310e+001  1.1938e+001  1.2566e+001  1.3195e+001  1.3823e+001...
    1.4451e+001  1.5080e+001  1.5708e+001  1.6336e+001  1.6965e+001  1.7593e+001  1.8221e+001  1.8850e+001  1.9478e+001  2.0106e+001  2.0735e+001...
    2.1363e+001  2.1991e+001  2.2619e+001  2.3248e+001  2.3876e+001  2.4504e+001  2.5133e+001  2.5761e+001  2.6389e+001  2.7018e+001  2.7646e+001...
    2.8274e+001  2.8903e+001  2.9531e+001  3.0159e+001  3.0788e+001  3.1416e+001  3.2044e+001  3.2673e+001  3.3301e+001  3.3929e+001  3.4558e+001...
    3.5186e+001  3.5814e+001  3.6442e+001  3.7071e+001  3.7699e+001  3.8327e+001  3.8956e+001  3.9584e+001  4.0212e+001  4.0841e+001  4.1469e+001...
    4.2097e+001  4.2726e+001  4.3354e+001  4.3982e+001  4.4611e+001];  % [rad/s]
m2=[ -2.4555e+001 -1.8863e+001 -1.5273e+001 -1.3022e+001 -1.1368e+001 -1.0106e+001 -9.0412e+000 -8.1258e+000 -7.5362e+000 -6.8305e+000 -6.3075e+000...
    -5.8365e+000 -5.3956e+000 -4.9887e+000 -4.6133e+000 -4.2957e+000 -3.9825e+000 -3.7159e+000 -3.5035e+000 -3.2543e+000 -3.0982e+000 -2.8406e+000...
    -2.5244e+000 -2.1561e+000 -1.9043e+000 -1.5786e+000 -1.3299e+000 -1.0794e+000 -1.0280e+000 -1.0049e+000 -1.0973e+000 -1.1484e+000 -1.3472e+000...
    -1.6581e+000 -1.7930e+000 -2.1920e+000 -2.0484e+000 -2.6549e+000 -2.3225e+000 -2.4648e+000 -2.8349e+000 -2.5328e+000 -2.5987e+000 -2.5688e+000...
    -3.6216e+000 -3.6520e+000 -3.9749e+000 -3.6494e+000 -3.8424e+000 -3.3791e+000 -3.0330e+000 -3.3493e+000 -3.0616e+000 -3.6999e+000 -3.2204e+000...
    -3.6160e+000 -3.4716e+000 -3.2962e+000 -3.5312e+000 -3.2992e+000 -3.5701e+000 -3.3806e+000 -3.7542e+000 -3.7192e+000 -3.5343e+000 -3.8899e+000...
    -3.9252e+000 -4.1289e+000 -4.7280e+000 -4.2376e+000 -4.8412e+000];  % [dB]
ph=[7.9373e+001  7.3436e+001  6.8009e+001  6.2013e+001  5.6244e+001  5.0607e+001  4.6156e+001  4.1887e+001  3.7741e+001  3.3938e+001  3.0157e+001...
    2.6500e+001  2.2878e+001  1.9852e+001  1.6389e+001  1.2625e+001  9.3336e+000  6.4858e+000  3.4874e+000  8.1324e-001 -2.2675e+000 -5.3157e+000...
    -8.9050e+000 -1.3065e+001 -1.5734e+001 -1.9750e+001 -2.2638e+001 -2.6510e+001 -2.9596e+001 -3.2418e+001 -3.5896e+001 -4.0109e+001 -4.3411e+001...
    -4.7108e+001 -5.0974e+001 -5.4830e+001 -5.9614e+001 -6.3986e+001 -6.6701e+001 -7.0808e+001 -7.3210e+001 -8.2986e+001 -8.3098e+001 -8.6777e+001...
    -9.1118e+001 -9.3465e+001 -9.4334e+001 -9.8264e+001 -9.8835e+001 -9.9147e+001 -1.0102e+002 -1.0385e+002 -1.0682e+002 -1.0833e+002 -1.0766e+002...
    -1.1003e+002 -1.0825e+002 -1.1121e+002 -1.0969e+002 -1.1306e+002 -1.1880e+002 -1.1594e+002 -1.2291e+002 -1.2404e+002 -1.2812e+002 -1.2867e+002...
    -1.3082e+002 -1.3201e+002 -1.3627e+002 -1.3741e+002 -1.4003e+002];

zfr = m2.*exp(i*ph*pi/180);
Ts = 0;
data = idfrd(zfr,f,Ts);

A = [1 1 1 1 1 -1; 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0];
B = [1; 0; 0; 0; 0; 0];
C = [1 1 1 1 1 0];
D = [0];

m = idss(A,B,C,D,'Ts',0);

m.As = [NaN NaN NaN NaN NaN -1; 1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0; 0 0 0 1 0 0; 0 0 0 0 1 0]
m.Bs = [1; 0; 0; 0; 0; 0];
m.Cs = [NaN NaN NaN NaN NaN 0];
m.Ds = [0];

res = pem(data,m);

A=res.a;
B=res.b;
C=res.c;
D=res.d;

[num1, den1] = ss2tf(A, B, C, D);

[m1,ph1,f1] = bode(num1,den1);

tf(num1, den1)

figure(1)
plot(f/(2*pi),m2,'o','LineWidth',2)
hold on
grid on
plot(f1/(2*pi),20*log10(m1(:)),'r','LineWidth',2)
axis([0 7.5 -26 50])
xlabel('Frequency [2*pi*f]','VerticalAlignment','top','fontsize',17)
ylabel('Magnitiude [dB]','VerticalAlignment','bottom','fontsize',17)
legend('Measured data','New',1)

figure(2)
plot(f/(2*pi),ph(:),'o','LineWidth',2)
hold on
plot(f1/(2*pi),ph1(:),'r','LineWidth',2)
grid on
axis([0 7.5 -150 100])
xlabel('Frequency [2*pi*f]','VerticalAlignment','top','fontsize',17)
ylabel('Phase [deg]','VerticalAlignment','bottom','fontsize',17)
legend('Measured data','New',1)

[Kamil]

Any solutions???? Anybody can help me?