EE342.01: MATLAB M-FILE FOR SIMULATING DISCRETE-TIME DC MOTOR CONTROL

%
%  Filename:  example7.m
%
%  Description:  M-file to demonstrate discrete-time control of a DC motor.
%

clear;                                  % clear matlab memory

tmax = 0.20;                            % amount of time to run simulation
T = 0.001;                              % sample interval

Kp = 700;                               % proportional controller gain
Kv = 0.4;                               % derivative controller gain
Ki = 1000;                               % integral controller gain
wd = 1.0;                               % desired angular velocity (krpm)
eint = 0;                               % initialize error integral to zero

KT = 4.07;				% motor torque constant
L = 0.00156;                            % motor inductance
J = 2.0;                                % motor rotor & load inertia
R = 1.89;                               % motor resistance
B = 13.0;                               % motor/gear/bearing friction
Kb = 3.01;                              % motor back emf constant

wprev1 = 0;                             % w[1] = wdot(0)*T + w(0)
wprev2 = 0;                             % w[0] = w(0) = 0

eprev1 = wd - wprev1;                   % initialize errors
eprev2 = wd - wprev2;

nTvec = [0 T]; wvec = [wprev2 wprev1];  % store ICs for plotting

for n=2:tmax/T,

  eint = eint + eprev2*T;               % sum (integrate) e
  eder = (eprev1 - eprev2)/T;           % differentiate e using Euler approx.
  
  vaprev2 = Kp*eprev2 + Kv*eder + Ki*eint;  % compute input voltage

  w = (2-(R*J+B*L)*T/(L*J))*wprev1 + ...
      ((R*J+B*L)*T/(L*J) - 1 - (B*R+KT*Kb)*T*T/(L*J))*wprev2 + ...
      KT*T*T/(L*J)*vaprev2;             % DC motor difference eqn

  wprev2 = wprev1;                      % store values for next iteration
  wprev1 = w;
  
  eprev2 = eprev1;
  eprev1 = wd - wprev1;
  
  nTvec(n+1) = n*T;                     % store nT,w into vectors for plotting
  wvec(n+1) = w;
  
end

plot(nTvec, wvec, 'o'); grid;
xlabel('time (sec)'); ylabel('angular velocity (krpm)');
ttle = ['Motor Step Response for Kp=', num2str(Kp), ', Kv=', num2str(Kv), ...
        ', Ki=', num2str(Ki)];
title(ttle);