| cholonco |
23-04-2013 04:28 PM |
Thiết kế hệ thống dùng matlab
Thiết kế hệ thống dùng matlab
Phương pháp sớm pha dùng giản đồ Bode, hàm leadbode
Mã:
function[gc]=leadbode(g,K,PM,GM)
%This function using the Bode diagram make the phase lead compensation,
%input is tf, Kp(or Kv, Ka) (if have no K, Kc=1),phase margin and gain
%margin which expect.
%-------------------------------------------------------|%
|%
%-------------------------------------------------------|%
if (nargin==3)
GM=PM;
PM=K;
Kc=1;
elseif (nargin==4)
Kc=Kc_error(g,K);
end
%Now, find margin
g1=Kc*g;
[gm1,Pm,wcg1,wcp1]=margin(g1);
clear gm1 wcp1 wcg1
%Step 4 in course, loop at here...
%Phase angle adder
%Step 5
theta=2;
max_ang=30;
while(theta
phimax=(PM-Pm+theta)*pi/180;
if (phimax>89*pi/180)
Result=sprintf('Bai toan vo nghiem')
gc=nan;
break;
end
%Step 6
anpha=(1+sin(phimax))/(1-sin(phimax));
%Step 7
[mag phase w]=bode(g1);
clear phase
mag=mag(:);
w=w(:);
%Find freq value which satisfies relation: Abs(G1(jw'c))=1/sqrt(anpha)
temp=1/sqrt(anpha);
t=find((mag>temp-7)&(mag
[vltemp,degree]=min(abs(mag(t)-temp));
clear vltemp mag
wc=w(min(t)+degree-1);
clear w degree t
wc=wc-0.5:0.001:0.5+wc;
mag=abs(freqresp(g1,wc));
mag=mag(:);
[value,degree]=min(abs(mag-temp));
clear value
wc=wc(degree);
T=1/wc*temp;
gc=Kc*tf([anpha*T 1],[T 1]);
%Step 9
gkq=g*gc;
[gmc pmc wgc wpc]=margin(gkq);
gmc=20*log10(gmc);
if ((gmc>=GM)&(pmc>=PM))
clf;
hold on
bode(gc,g);
bode(gkq,'r')
legend('Phase lead compensation','Org Tranfer Function','Compensated Tranfer Function');
margin(g*gc);
pause();
clf;
hold on;
step(feedback(gkq,1),feedback(g,1));
grid on
legend('Compensated Tranfer Function','Org Tranfer Function');
pause();
close();
break
else
clear gc gmc pmc wgc wpc gkq
theta=theta+1;
if (theta==max_ang)
Result=sprintf('Bai toan vo nghiem')
gc=nan;
end
end
end
%End
Phương pháp trễ pha dùng giản đồ Bode, hàm lagbode
Mã:
function[gc]=lagbode(g,K,PM,GM)
%This function make the phase lag compensation, input is tf,Kp(or Kv,Ka)
%phase margin and gain margin which expect.(Kc=1) in this func if have no K.
%gc=lagbode(g,K,PM,GM)
%-------------------------------------------------------|%
|%
%-------------------------------------------------------|%
if (nargin==3)
Kc=1;
%This code to find Kc
elseif (nargin==4)
Kc=Kc_error(g,K);
end
g1=Kc*g;
[gm1,Pm,wcg1,wcp1]=margin(g1);
clear gm1 wcp1 wcg1
%Phase angle adder
theta=2;
max_ang=30;
%Step 3
while(theta
phi=(-180+PM+theta);
%Find freq value which satisfies relation: Phi1((w'c))=-180+PM+theta
[mag,phase,w]=bode(g1);
clear mag
phase=phase(:);
w=w(:);
t=find((phase>phi-10)&(phase
[vltemp,degree]=min(abs(phase(t)-phi));
if (degree>=1)
clear vltemp phase
wc=w(min(t)+degree-1);
clear w degree t
w=abs(wc-0.5):0.001:(0.5+wc);
w=w(:);
[mag,phase]=bode(g1,w);
phase=phase(:);
[value,degree]=min(abs(phase-phi));
wc=w(degree);
clear value w
%Find anpha
anpha=1/abs(freqresp(g1,wc));
%Caculate time constant T
T=ceil(12/(wc*anpha));
clear wc
%Phase lag compensation
gc=Kc*tf([anpha*T 1],[T 1]);
%Step 7
gkq=g*gc;
[gmc,pmc,wgc,wpc]=margin(gkq);
gmc=20*log10(gmc);
if ((gmc>GM)&(pmc>PM))
clf;
hold on
bode(gc);
bode(g);
bode(gkq,'r')
legend('Phase lag compensation','Org Tranfer Function','Compensated Tranfer Function');
margin(gkq);
pause();
clf;
if(g.iodelay==0)
clf;
hold on;
step(feedback(gkq,1),feedback(g,1));
grid on
legend('Compensated Tranfer Function','Org Tranfer Function');
pause();
close();
end
break
else
theta=theta+1;
if (theta==max_ang)
Result=sprintf('Bai toan vo nghiem')
gc=nan;
end
end
else
clear gc gmc pmc wgc wpc gkq
theta=theta+1;
if (theta==max_ang)
gc=nan;
Result=sprintf('Bai toan vo nghiem')
end
end
end
%End
Bạn cần load hàm Kc_error để có thể chạy được các hàm trên.
Mã:
function[Kc]=Kc_error(g,k)
%Kc=Kc_error(g,k)
%Assume that:
%If tf have no ideal integral, k=Kp
%If tf have 1 ideal integral, k=Kv
%If tf have 2 ideal integral, k=Ka
%-------------------------------------------------------|%
%|(C)2005 Bui Trung Hieu |%
%|Website: www.khvt.com |%
%|Email: hieu@khvt.com |%
%-------------------------------------------------------|%
if (nargin==2)
if(length(k)==1)
sw=length(find(pole(g)==0));
if (sw==0)
%This code to find Kc from Kp
a=dcgain(g);
Kc=double(k/a);
elseif (sw==1)
%This code to find Kc from Kv
s=tf('s');
a=dcgain(minreal(s*g));
Kc=double(k/a);
elseif (sw==2)
%This code to find Kc from Ka
s=tf('s');
a=dcgain(minreal(s^2*g));
Kc=double(k/a);
end
else
Kc=[];
error('???Ban nhap sai so chieu cua k!!!');
end
else
Kc=[];
error('???Xem help de co them thong tin!!!');
end
Nguon: Suu tam
|
