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

其他行业

开发平台:

Matlab

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