clear

clc

//Example 13.8

disp('Example 13.8')

s = %s;

num = 4;

den = (5*s+1);

delay=1;

w = 0.001:0.002:10*%pi;

LF = "loglog" // Warning: Change this as necessary

Gv=2;Gm=0.25;

Ku=4.25;Pu=2*%pi/1.69;

//Ziegler Nichols

Kc1=0.6*Ku;taui1=Pu/2;tauD1=Pu/8;

//Tyreus Luyben

Kc2=0.45*Ku;taui2=Pu*2.2;tauD2=Pu/6.3;

mprintf(' Kc tauI tauD')

mprintf('\nZN %f %f %f',Kc1,taui1,tauD1)

mprintf('\nTL %f %f %f',Kc2,taui2,tauD2)

Gc_ZN=Kc1*(1+1/taui1/s+s*tauD1/(0.1*s*tauD1+1));

Gc_TL=Kc2*(1+1/taui2/s+s*tauD2/(0.1*s*tauD2+1)); //Filtered Controllers with filter constant as 0.1

G1 = num/den*Gc_ZN*Gm*Gv;

G1m = horner(G1,%i*w); //G1m denotes magnitude

Abs_G1m=abs(G1m)

G1p = phasemag(G1m)-delay*w*180/%pi; //G1p denotes phase

G2 = num/den*Gc_TL*Gm*Gv;

G2m = horner(G2,%i*w); //G1m denotes magnitude

Abs_G2m=abs(G2m);

G2p = phasemag(G2m)-delay*w*180/%pi; //G1p denotes phase

xset('window',0);

subplot(2,1,1)

plot2d(w,Abs_G1m,logflag="ll",style = 1,rect=[0.01,0.1,10,1000]);

plot2d(w,Abs_G2m,logflag="ll",style = 2,rect=[0.01,0.1,10,1000]);

legend("Ziegler-Nichols","Tyreus Luyben")

xgrid();

xtitle('Loglog','','Magnitude');

subplot(2,1,2)

plot2d1(w,G1p,logflag="ln",style = 1,rect=[0.01,-300,10,-50]);

plot2d1(w,G2p,logflag="ln",style = 2,rect=[0.01,-300,10,-50]);

// legend("Ziegler-Nichols","Tyreus Luyben")

xgrid();

xtitle('w(rad/sec)','','Phase(deg)');

//G_ZN=syslin('c',G1);

//G_ZN=iodelay(G_ZN,delay);

//G_TS=syslin('c',G2);

//G_TS=iodelay(G_TS,delay);

//scf();nyquist(G_TS,%f)

//[gm_ZN,fr_ZN]=g_margin(G_ZN);[gm_TS,fr_TS]=g_margin(G_TS);

//[pm_ZN,fr_ZN_p]=p_margin(G_ZN);[pm_ZN,fr_ZN_p]=p_margin(G_TS);

//g_maring and p_margin do not support iodelay toolbox, hence we

//cannot use these so we try a workaround

//We can find w for which magnitude(AR) is 1 and

//and calculate phase corresponding to it which gives us phase margin

//Also we can find crossover frequency and thus find Gain Margin

indices1=find(abs(Abs_G1m-1)

Two rather esoteric methods of frequency component estimation for complex spectra are given below.  They describe the use of the harmonic product spectrum (Rabiner and Gold), and the heterodyne filter.

https://books.google.com/books/about/Theory_and_application_of_digital_signal.html?id=iAxTAAAAMAAJ

https://books.google.com/books?id=wIZt9rZ27i4C&pg=PA5&lpg=PA5&dq=heterodyne+filter+analysis&source=bl&ots=A9-eVYpoWK&sig=gITMjs63wXbspgCufkW-lZ1IguM&hl=en&sa=X&ved=0ahUKEwjJstW7jMjLAhUCwWMKHVyEBfMQ6AEILjAC#v=onepage&q=heterodyne%20filter%20analysis&f=false

Hello Hadi,

I do not have any software code for your requirement but there is one method which I thought might be easier for you to understand and implement.

Kindly refer book on Power Electronics by Theodore Wildi, 3rd Edn. Look for chapter on Harmonic Analysis. I think it would be best suited to you.

Best wishes!

-Prasanna Waichal

Dear Hadi,

Hope that this matlab file will give you the answer for your queries...

Regards