home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 10.ddi / CONTROL.DI$ / LSIM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.7 KB  |  118 lines
  1. function  [yout,x] = lsim(a, b, c, d, u, t, x0)
  2. %LSIM    Simulation of continuous-time linear systems to arbitrary inputs.
  3. %    LSIM(A,B,C,D,U,T) plots the time response of the linear system:
  4. %            .
  5. %            x = Ax + Bu
  6. %            y = Cx + Du
  7. %    to the input time history U. Matrix U must have as many columns as
  8. %    there are inputs, U.  Each row of U corresponds to a new time 
  9. %    point, and U must have LENGTH(T) rows.  The time vector T must be
  10. %    regularly spaced.  LSIM(A,B,C,D,U,T,X0) can be used if initial 
  11. %    conditions exist.
  12. %
  13. %    LSIM(NUM,DEN,U,T) plots the time response of the polynomial 
  14. %    transfer function  G(s) = NUM(s)/DEN(s)  where NUM and DEN contain
  15. %    the polynomial coefficients in descending powers of s.  When 
  16. %    invoked with left hand arguments,
  17. %        [Y,X] = LSIM(A,B,C,D,U,T)
  18. %        [Y,X] = LSIM(NUM,DEN,U,T)
  19. %    returns the output and state time history in the matrices Y and X.
  20. %    No plot is drawn on the screen.  Y has as many columns as there 
  21. %    are outputs, y, and with LENGTH(T) rows.  X has as many columns 
  22. %    as there are states.
  23. %
  24. %    See also: STEP,IMPULSE,INITIAL and DLSIM.
  25.  
  26. %    LSIM normally linearly interpolates the input (using a first order hold)
  27. %    which is more accurate for continuous inputs. For discrete inputs such 
  28. %    as square waves LSIM tries to detect these and uses a more accurate 
  29. %    zero-order hold method. LSIM can be confused and for accurate results
  30. %    a small time interval should be used.
  31.  
  32. %    J.N. Little 4-21-85
  33. %       Revised A.C.W.Grace 8-27-89 (added first order hold)
  34. %                        1-21-91 (test to see whether to use foh or zoh)
  35. %    Copyright (c) 1986-93 by the MathWorks, Inc.
  36.  
  37. error(nargchk(4,7,nargin));
  38.  
  39. if (nargin==4),            % transfer function description 
  40.     [num,den] = tfchk(a,b);
  41.     u = c;
  42.     t = d;
  43.     % Convert to state space
  44.     [a,b,c,d] = tf2ss(num,den);
  45.  
  46. elseif (nargin==5),
  47.      error('Wrong number of input arguments.');
  48. else
  49.     error(abcdchk(a,b,c,d));
  50. end
  51.  
  52. [ns,n] = size(a);
  53. if (nargin==6)|(nargin==4),
  54.     x0 = zeros(ns,1);
  55. end
  56.  
  57. [p,m]=size(d);
  58. if p*m==0, x=[]; if nargout~=0, yout=[]; end; return; end
  59.  
  60. if m==1, u=u(:); end, [nu,mu]=size(u);
  61. t=t(:); nt = length(t);
  62.  
  63. % Make sure u has the right number of columns and rows.
  64. if mu~=m, error('U must have the same number of columns as inputs.'); end
  65. if nu~=nt, error('U must have the same number of rows as the length of T.'); end
  66.     
  67. TS=t(2)-t(1);
  68. % First Order Hold Approximation 
  69.  
  70. % Try to find out whether the input is step-like or continuous
  71.  
  72. % Differentiate u
  73. udiff = u(2:nu,:)-u(1:nu-1,:);
  74.  
  75. % Find all transition points of the input i.e., all changing inputs
  76. uf = find(udiff ~= 0);
  77.  
  78. % Find out if changes in the input occur after a number of constant inputs.
  79. % If they do then assume that it is a step-like input and use a zero 
  80. % order hold method. Otherwise assume that the input is continuous and 
  81. % linearly interpolate using first order hold method.
  82. usteps = find(diff(uf) == 1);
  83. foh = (length(usteps) > 0);
  84.  
  85. if foh 
  86.     % Use first order hold approximation 
  87.  
  88. % For first order hold approximation first add m integrators in series
  89.     [a,b,c,d]=series(zeros(m),eye(m),eye(m),zeros(m),a,b,c,d);
  90.  
  91. % For first order hold add (z-1)/TS in series
  92. % This is equivalent to differentiating u.
  93. % Transfer first sample to initial conditions.
  94.  
  95.     x0=[zeros(m,1);x0(:)]+(b*u(1,:).');
  96.     u(1:nu-1,:) = udiff/TS;
  97.     u(nu,:)=u(nu-1,:);
  98. end
  99.  
  100. % Get equivalent zero order hold discrete system
  101. [A,B] = c2d(a,b,t(2)-t(1));
  102.  
  103.  
  104. x = ltitr(A,B,u,x0);
  105. y = x * c.' + u * d.';
  106.  
  107. if foh
  108.     % Remove the integrator state
  109.     x=x(:,1+m:n+m);
  110. end
  111.  
  112. if nargout==0        % If no output arguments, plot graph
  113.     plot(t,y), xlabel('Time (secs)'), ylabel('Amplitude')
  114.     return % Suppress output
  115. end
  116. yout = y; 
  117.  
  118.