home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / DISKDEMO.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.5 KB  |  184 lines
  1. %DISKDEMO.M Demonstration design of harddisk digital controller.
  2. echo on
  3.  
  4. % This file demonstrates MATLAB's ability for classical digital control
  5. % system design by going through the design of a computer HARDDISK 
  6. % read/write head position controller.
  7.  
  8. echo off
  9. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  10. echo on
  11.  
  12. pause % Press any key to continue ...
  13.  
  14. % Using Newton's law, the simplest model for the read/write head has the 
  15. % following differential equation:
  16. %
  17. %  I*theta_ddot + C*theta_dot + K*theta = Ki * i
  18. %
  19. %  where I is the inertia of the head assembly
  20. %        C is the viscous damping coefficient of the bearings
  21. %        K is the return spring constant
  22. %        Ki is the motor torque constant
  23. %        theta_ddot, theta_dot, and theta are the angular acceleration,
  24. %          angular velocity and position of the head
  25. %    i is the input current
  26. %
  27.  
  28. pause % Press any key to continue ...
  29.  
  30. % Taking the laplace transform, the transfer function is:
  31. %                Ki
  32. %  H(s) = ----------------
  33. %         I s^2 + C s + K
  34. %
  35. % Using the values I=.01 Kg m^2, C=.004 Nm/(rad/sed), K=10 Nm/rad,  and
  36. % Ki=.05 Nm/rad form the transfer function description of this system
  37.  
  38. I = .01; C = 0.004; K = 10; Ki = .05;
  39. NUM = [Ki];
  40. DEN = [I C K];
  41. printsys(NUM,DEN,'s');
  42.  
  43. pause % Press any key to continue ...
  44.  
  45. % Our task is to design a digital controller that can be used to provide 
  46. % accurate positioning of the read/write head.  We will do the design in the
  47. % digital domain.  
  48.  
  49. % First we must discretize our plant since it is continuous.  Since our
  50. % plant will have a digital-to-analog-converter (with a zero-order hold)
  51. % connected to its input, use the 'zoh' discretization method
  52. % of the function C2DM. Use sample time Ts = 0.005  (5 ms)
  53.  
  54. Ts = 0.005;
  55. w = logspace(0,3);
  56. [mag,phase] = bode(NUM,DEN,w);
  57. [num,den] = c2dm(NUM,DEN,Ts,'zoh');
  58. [mzoh,pzoh,wzoh] = dbode(num,den,Ts);
  59.  
  60. % Now plot the results as a comparison.  Press any key after the plot ...
  61. echo off
  62. subplot(211)
  63. semilogx(w,20*log10(mag),wzoh,20*log10(mzoh))
  64. xlabel('Frequency (rad/sec)'), ylabel('Gain db')
  65. title('c2d comparison plot')
  66.  
  67. subplot(212)
  68. semilogx(w,phase,wzoh,pzoh)
  69. xlabel('Frequency (rad/sec)'), ylabel('Phase deg')
  70. pause % Press any key to continue ...
  71. echo on
  72.  
  73. pause % Press any key to continue ...
  74.  
  75. % Now analyze the discrete system.
  76. disp('Discrete system')
  77. printsys(num,den,'z')
  78.  
  79. % Plot step response
  80. subplot(111)
  81. dstep(num,den); pause % Press any key after the plot ...
  82.  
  83. % The system oscillates quite a bit.  This is probably due to very light
  84. % damping.  We can check this by computing the open loop eigenvalues.
  85.  
  86. disp('Open loop discrete eigenvalues'), ddamp(den,Ts); 
  87. zgrid('new'), pzmap(1,den); pause % Press any key after the plot ...
  88.  
  89. % Note that the poles are very lightly damped and near the unit circle.
  90. % We need to design a compensator that increases the damping of this system.
  91.  
  92. % Let's try to design a compensator.  The simplest compensator is a simple gain.
  93. rlocus(num,den); hold off; pause % Press any key after the plot ...
  94.  
  95. % As shown in the root locus, the poles quickly leave the unit circle and go
  96. % unstable.  We need to introduce some lead or a compensator with some zeros.
  97. % Try the compensator:        K(z + a)
  98. %                     D(z) = --------  where a < b
  99. %                             (z + b)
  100.  
  101. pause % Press any key to continue ...
  102.  
  103. % Form compensator and connect in series with our plant
  104. % Use a = -.85 and b = 0.
  105. [numc,denc] = zp2tf([.85 ]',[0]',1);
  106. [num2,den2] = series(numc,denc,num,den);
  107.  
  108. % Lets see how the frequency response has changed.
  109. [mag,phase,w] = dbode(num,den,1);    % Use normalized frequency
  110. [mag2,phase2] = dbode(num2,den2,1,w);
  111.  
  112. % Now plot a comparison plot.  Press any key after the plot ...
  113. echo off
  114. subplot(211), semilogx(w,20*log10(mag),w,20*log10(mag2))
  115. xlabel('Frequency (rad/sec)'), ylabel('Gain dB')
  116. subplot(212), semilogx(w,phase,w,phase2)
  117. xlabel('Frequency (rad/sec)'), ylabel('Phase deg')
  118. % Plot -180 degree line
  119. hold on; semilogx([min(w(:)),max(w(:))],[-180,-180],'w--'); hold off;
  120. pause % Press any key to continue ...
  121. echo on
  122.  
  123. % So our compensator has indeed added lead.
  124.  
  125. % Now let's try the root locus again with our compensator
  126. subplot(111)
  127. zgrid('new'), rlocus(num2,den2); hold off; pause % Press any key after plot ...
  128.  
  129. % This time the poles stay within the unit circle for some time.
  130. % Now its your turn, Using RLOCFIND chose the poles with greatest damping
  131. % ratio.  (The lines drawn by ZGRID show the damping ratios from z=0 to 1
  132. % in steps of .1)
  133.  
  134. pause % Press any key and then choose a point on the plot
  135. [k,poles] = rlocfind(num2,den2);
  136.  
  137. disp(['You chose gain: ',num2str(k)]), ddamp(poles,Ts);
  138.  
  139. % Let's form the closed loop system so that we can analyze the design.
  140. [numc,denc] = feedback(num2,den2,k,1);
  141.  
  142. % These eigenvalues should match the ones you chose.
  143. disp('Closed loop eigenvalues'), ddamp(denc,Ts);
  144.  
  145. pause % Press any key to continue ...
  146.  
  147. % Closed loop time response
  148. dstep(numc,denc); pause % Press any key after the plot ...
  149.  
  150. % So the response looks pretty good and settles in about 14 samples
  151. % which is 14*Ts secs.
  152.  
  153. disp(['Our disc drive will have a seek time > ',num2str(14*Ts),' seconds.'])
  154.  
  155. pause % Press any key to continue ...
  156.  
  157. % Let's now look at the robustness of our design.  The most common classical
  158. % robustness criteria is the gain and phase margin.  The criteria is determined
  159. % by forming a unity feedback system, calculating the Bode response and looking
  160. % for the phase and gain crossovers.  MATLAB contains a function MARGIN that
  161. % determines the phase and gain margin given the Bode response.
  162.  
  163. % Form unity feedback system by connecting our design system with the gain
  164. % we chose.  Leave the loop open so we can compute the open loop Bode response.
  165. [num2,den2] = series(num2,den2,k,1);
  166.  
  167. % Compute Bode response and margins
  168. [mag,phase,w] = dbode(num2,den2,Ts);
  169. [Gm,Pm,Wcg,Wcp] = margin(mag,phase,w); 
  170.  
  171. % Plot Bode plot with margins
  172. margin(mag,phase,w); pause % Press any key after the plot ...
  173.  
  174. % Gain margin db, @ frequency, Phase margin, @ frequency
  175. Margins = [20*log10(Gm),Wcg,Pm,Wcp]
  176.  
  177. % Our design is robust and can tolerate a 10 db gain increase and a 40 degree
  178. % phase lag without going unstable.  By continuing this design process we may
  179. % be able to find a compensator that will stabilize the open loop system and 
  180. % allow us to reduce the seek time (more damping would allow us to reduce the
  181. % settling time in the step response).
  182.  
  183. echo off
  184.