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

其他行业

开发平台:

Matlab

  1. function [numopen, denopen, denclsd] = phlag(num,den,zeta)
  2. % Hadi Saadat, 1998
  3. discr=[
  4. '                                                                         '
  5. '  The function phlag(num, den, zeta) is used for the design of a phase-  '
  6. '  lag compensator using the approximated method k =~ K0/Kc.  num & den   '
  7. '  are row vectors of polynomial coefficients of the uncompensated open-  '
  8. '  loop plant transfer function.  zeta is the desired damping ratio. The  '
  9. '  user is prompted to enter gain K required for the steady-state error   '
  10. '  specification & the magnitude of the compensator zero, Z0. [Z0 <<硈1砞.'
  11. '  The function returns the open-loop and the closed-loop numerator and   '
  12. '  denominator of the compensated system transfer function.               '
  13. '                                                                         '];
  14. %disp(discr);
  15. K=input('Enter gain K required for the steady-state error specification -> ');
  16. Z0=input('Enter magnitude of the compensator zero -> ');
  17. clc
  18. %beta=atan2(imag(s1),real(s1));
  19. %sm=abs(s1);
  20. %[mag,phase]=ghs(num,den,s1);   % Returns the mag. and phase of G(s1)H(s1)
  21. c1=0:.1:10; c2=11:1:100; c3=110:10:1000;
  22. k=[c1 c2 c3]; kk=k';
  23. j=0;
  24. r=rlocus(num,den,k);
  25.  ri=imag(r);  risum=sum(ri);
  26. costhta=cos([pi-angle(r)]);
  27. j= find (risum > .5);
  28. if isempty(j)==1;
  29.        disp('Warning: zeta line does not cross root-locus')
  30.        return,else, end
  31.  l=length(costhta);
  32.  for i=1:l for c=1:j
  33.      if costhta(i,c) ==-1 costhta(i,c)=1; end,end,end
  34. cost=(costhta(:,j));
  35. l=length(cost);
  36. if cost(1) > zeta
  37.         i = find(cost < zeta);
  38. else
  39.         i = find(cost > zeta);
  40. end
  41. if length(i) == 0
  42.         disp('Warning: zeta line does not cross root-locus'),return
  43. else
  44. while i(1)>l i(1)=i(1)-l; end
  45.         i2 = i(1);
  46.         i1 = i2 - 1;
  47.         k1 = min(kk(i1,:));
  48.         k2 = min(kk(i2,:)) ;
  49.         m1 = min(cost(i1,:));
  50.         m2 = min(cost(i2,:));
  51.         K0 = k1 + (m1-zeta)*(k2-k1)/(m1-m2);
  52. end
  53. K0=round(100*K0)/100;
  54. clc
  55. fprintf('   Gain for the desired closed-loop pole      K0 = %gnn',K0)
  56. fprintf('   Gain for the desired steady-state response K  = %gnn',K)
  57. Kc=K0/K;
  58. P0=Kc*Z0;
  59. KKc=K*Kc;
  60. fprintf('   Gc(0) = %g',1),fprintf(',      Gc = %g',Kc),fprintf('(s + %g',Z0),
  61. fprintf(')/(s + %g',P0), fprintf(')nn')
  62. % the following statements will form the characteristic Equation
  63. % of the compensated system.
  64. m=length(num); n=length(den);
  65. if n > m
  66. o=zeros(1,n-m); mk=[o,1]; num1=conv(num,mk);
  67. else, num1=num, end
  68. numgc=[KKc,KKc*Z0]; numopen=conv(numgc,num1);
  69. dengc=[1,P0];       denopen=conv(dengc,den);
  70. denclsd=denopen+numopen;
  71. %fprintf('Row vectors of polynomial coefficients of the compensated system:nn')
  72. %fprintf('Open-loop num. '),disp(numopen)
  73. %fprintf('Open-loop den. '),disp(denopen)
  74. %fprintf('Closed-loop den'),disp(denclsd)
  75. fprintf('Compensated open-loop ')
  76. GH = tf(numopen, denopen)
  77. fprintf('Compensated closed-loop ')
  78. T = tf(numopen, denclsd)
  79. discr2=[
  80. 'Roots of the compensated characteristic equation:       '];
  81. disp(discr2)
  82. r=roots(denclsd);
  83. disp(r)
  84. rreal=real(r);
  85.   for l=1:n
  86.   if rreal(l) >=0
  87.   fprintf('   Root on the RHP, system is unstable.n')
  88.   fprintf('   Change Z0 or K & repeat.nn')
  89.   else,end
  90.   end