Frqlag.m
上传用户:eighthdate
上传日期:2014-05-24
资源大小:270k
文件大小:5k
源码类别:

其他行业

开发平台:

Matlab

  1. function [numopen, denopen, denclsd] = frqlead(num, den)
  2. % Hadi Saadat, 1998
  3. %discr=[
  4. %'                                                                           '
  5. %'  The function [numopen,denopen,denclsd]=frqlag(num, den) is used for the  '
  6. %'  frequency response design of a phase-lag compensator.  num & den are     '
  7. %'  row vectors of polynomial coefficients of the uncompensated open-loop    '
  8. %'  plant transfer function. Design is based on the phase margin criterion.  '
  9. %'  The user is prompted to enter the desired Phase Margin & the controller  '
  10. %'  DC gain Gc(0) = Kc*Z0/P0.  The program finds and displays a compensated  '
  11. %'  gain crossover frequency range for an stable controller.  The user is    '
  12. %'  then prompted to enter the gain crossover frequency in this range.  The  '
  13. %'  controller transfer function and the frequency-domain specifications     '
  14. %'  before and after compensation are found.  The function returns the open- '
  15. %'  loop and closed-loop numerator and denominators of the compensated system'
  16. %'  transfer function.                                                       '];
  17. %disp(discr);
  18. r=abs(roots(den));
  19. i=find(r>0); rp=r(i);
  20. rmx=max(rp); rmn=min(rp); wst=0.1*round(rmn); wf=20*round(rmx);;dw=wf/800;
  21. w=wst:dw:wf;
  22. clc
  23. a0=input('Enter the compensator DC Gain -> ');
  24. pm=input('Enter desired Phase Margin -> ');
  25. fprintf('n')
  26. [mag, phase] = bode(num, den, w);
  27. phase=180/pi*unwrap(phase*pi/180);
  28. if phase(1) > (-180+pm)
  29.      i = find(phase < (-180 + pm));
  30. else
  31.      i = find(phase > (-180 + pm));
  32. end
  33. if length(i)==0
  34.      disp('Phase does not cross (-180 + P.M.). ')
  35. return, else
  36.      i2=i(1); i1 = i2 -1;
  37.      if i1==0 wpm=w(i2); else
  38.      wa = w(i1); wb = w(i2); p1 = phase(i1); p2 = phase(i2);
  39.      wpm = wa + (-180+pm - p1)/(p2-p1)*(wb-wa);
  40.      end
  41. end
  42. i=find(mag<1/a0);
  43. if length(i)==0
  44.      disp('Gain does not cross 1/a0. ')
  45. return, else
  46.      i2=i(1); i1=i2 -1;
  47.      if i1==0 wga=w(i2); else
  48.      wa = w(i1); wb = w(i2); m1 = mag(i1); m2 = mag(i2);
  49.      wga =wa+(1/a0-m1)/(m2-m1)*(wb-wa);
  50.      end
  51. end
  52. w1=min(wga, wpm);
  53. [M,ph]=bode(num, den, w1);  % Returns the mag. and phase of G(w)H(w1)
  54. thta=-180 + pm - ph;
  55. thtar=thta*pi/180;
  56. [a1,b1]=frcntlr(num,den,w1,a0,pm);
  57. stab=0;
  58. wmn=w1/10; dw2=w1/10;
  59. while a1 < 0 | b1 < 0  & w1 > wmn
  60.   w1=w1-dw2;
  61.   [M,ph]=bode(num, den, w1);  % Returns the mag. and phase of G(w)H(w1)
  62.   thta=-180 + pm - ph;
  63.   thtar=thta*pi/180;
  64.   [a1,b1]=frcntlr(num,den,w1,a0,pm);
  65. end
  66. w1mx=w1;
  67. k=0;
  68. while a1 > 0 & b1 > 0  & w1 > wmn
  69.   k=k+1;
  70.   stab=stab+1;
  71.   w1=w1-dw2;
  72.   [M,ph]=bode(num, den, w1);  % Returns the mag. and phase of G(w)H(w1)
  73.   thta=-180 + pm - ph;
  74.   thtar=thta*pi/180;
  75.   [a1,b1]=frcntlr(num,den,w1,a0,pm);
  76.     if w1>wmn
  77.     wm=sqrt(a0/a1*1/b1); wdd = w1-wm;
  78.     wddm(k)=abs(wdd);  wk(k)=w1;
  79.     end
  80. end
  81.   [wmin,i]=min(wddm); wcal=wk(i);
  82. w1mn=w1;
  83. if stab==0
  84.   fprintf('Unstable controller. Reduce Phase Margin and repeat.n'),
  85.   return
  86.   else
  87.   fprintf('For a stable controller select a compensated gain crossovern')
  88.   fprintf('frequency wgc between %7.3g',w1mn),fprintf(' and %7.3gn',w1mx)
  89.     if wcal < wmn
  90.     fprintf('Suggested wgc is %7.3gn',wcal)
  91.     else,end
  92. end
  93. w1=input('Enter wgc -> ')
  94. clc
  95. fprintf('Uncompensated control system n')
  96. [Gm1, Pm1, wpc1, wgc1]=margin(mag, phase, w);
  97. fprintf('Gain Margin  = %7.3g',Gm1),fprintf('    Gain crossover w = %7.3gn',wgc1)
  98. fprintf('Phase Margin = %7.3g',Pm1),fprintf('   Phase crossover w = %7.3gn',wpc1)
  99. fprintf('n')
  100. [a1,b1]=frcntlr(num,den,w1,a0,pm);
  101. KKc=a1/b1; Z0=a0/a1; P0=1/b1;
  102. fprintf('Controller transfer function n')
  103. fprintf('   Gc(0) = %g',a0),fprintf(',     Gc = %g',KKc),fprintf('(s + %g',Z0),
  104. fprintf(')/(s + %g',P0), fprintf(')nn')
  105. % the following statements will form the characteristic Equation
  106. % of the compensated system.
  107. m=length(num); n=length(den);
  108. if n > m
  109. o=zeros(1,n-m); mk=[o,1]; num1=conv(num,mk);
  110. else, num1=num, end
  111. numgc=[KKc,KKc*Z0]; numopen=conv(numgc,num1);
  112. dengc=[1,P0];       denopen=conv(dengc,den);
  113. denclsd=denopen+numopen;
  114. %fprintf('Row vectors of polynomial coefficients of the compensated system:n')
  115. %fprintf('Open-loop num. '),disp(numopen)
  116. %fprintf('Open-loop den. '),disp(denopen)
  117. %fprintf('Closed-loop den'),disp(denclsd)
  118. fprintf('Compensated open-loop ')
  119. GH = tf(numopen, denopen)
  120. fprintf('Compensated closed-loop ')
  121. T = tf(numopen, denclsd)
  122. [magp,phasep]=bode(numopen,denopen,w);
  123. [Gm, Pm, wpc,wgc]=margin(magp,phasep,w);
  124. if Pm>360 Pm=Pm-360; else,end
  125. fprintf('Gain Margin  = %7.3g',Gm),fprintf('    Gain crossover w = %7.3gn',wgc)
  126. fprintf('Phase Margin = %7.3g',Pm),fprintf('   Phase crossover w = %7.3gn',wpc)
  127. fprintf('n')
  128. [M,ph]=bode(numopen, denclsd, w);
  129. frqspec(w,M)
  130. discr2=[
  131. 'Roots of the compensated characteristic equation:       '];
  132. disp(discr2)
  133. r=roots(denclsd);
  134. disp(r)
  135. rreal=real(r);
  136.    for l=1:n
  137.    if rreal(l) >=0
  138.   fprintf('   Root on the RHP, system is unstable.n')
  139.   else,end
  140.   end