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

其他行业

开发平台:

Matlab

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