home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / TRSP.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.3 KB  |  221 lines
  1. %function y = trsp(sys,input,tfinal,int,x0)
  2. %
  3. %  Calculates the time response of a SYSTEM matrix:
  4. %      SYS:    SYSTEM or CONSTANT matrix.
  5. %      INPUT:    VARYING input vector or constant.
  6. %      TFINAL:    final time value
  7. %              (optional, default = input final time).
  8. %      X0:    initial state
  9. %              (optional, default = 0).
  10. %      INT:    integration step
  11. %          (optional, default based on input & SYS. it
  12. %           is recommended that the user specify this,
  13. %           since the interval selected by TRSP is often
  14. %           quite small, and can make the calculation
  15. %           take a long time).
  16. %
  17. %   If INT is specified as 0 TRSP works out a suitable
  18. %   integration step. The output is generated at the selected
  19. %   integration step size.  To reduce the number of output points
  20. %   the user is referred to the function VDCMATE.
  21. %
  22. %   See also: DTRSP, SAMHLD, TUSTIN, VINTERP, and VDCMATE.
  23.  
  24. function y = trsp(sys,input,tfinal,int,x0)
  25.  
  26. if nargin == 0,
  27.     disp('usage: y = trsp(sys,input,tfinal,int,x0)')
  28.     return
  29.     end
  30.  
  31. [type,ny,nu,nx] = minfo(sys);
  32. if type == 'syst',
  33.     if nargin < 5,
  34.         x0 = zeros(nx,1);
  35.     else
  36.         [m,n] = size(x0);
  37.         if m == 1 & n == 1,
  38.              if x0 == 0,
  39.                  x0 = zeros(nx,1);
  40.                  [m,n] = size(x0);
  41.                  end
  42.              end
  43.         if m ~= nx | n ~= 1,
  44.             error('incorrect initial state vector dimension')
  45.             return
  46.                end
  47.         end
  48.     [a,b,c,d] = unpck(sys);
  49. elseif type ~= 'cons',
  50.     error('the system is not in packed or constant form')
  51.     return
  52.     end
  53.  
  54. %       check that the input vector is compatible and extract
  55. %       the data.
  56.  
  57. [typeu,nru,ncu,nupts] = minfo(input);
  58. if typeu == 'vary' & ncu == 1,
  59.     [udat,uptr,ut] = vunpck(input);
  60. elseif typeu == 'cons' & ncu == 1,
  61.     nupts = 1;
  62.     udat = input;
  63.     ut = 0;
  64. else
  65.     error('invalid input time vector')
  66.     return
  67.     end
  68.  
  69. if nru ~= nu, 
  70.     error('incompatible input vector dimension')
  71.     return
  72.     end
  73.  
  74. %       if an integration step size is supplied interpolate the
  75. %       input to this stepsize.  If one is not supplied then calculate
  76. %       a stepsize.
  77.  
  78. if nargin < 4,
  79.     int = 0;
  80.     end
  81. intstep = max(0,int);
  82. if intstep == 0,
  83.     if nupts > 1,
  84.         utnzdiff = diff(ut);
  85.         utnzdiff = utnzdiff(find(abs(utnzdiff) > eps));
  86.         if isempty(utnzdiff),
  87.             minstep = 1;
  88.         else
  89.             minstep = min(abs(utnzdiff));
  90.             end
  91.     else
  92.         minstep = 1;
  93.         end
  94.     if type == 'syst'
  95.         maxeig = max(abs(eig(a)));
  96.         intstep = min(minstep,0.1/maxeig);
  97.     else
  98.         intstep = minstep;
  99.         end
  100.     end
  101.  
  102. if intstep ~= int,
  103.     disp(['integration step size: ' num2str(intstep)])
  104.     end
  105.  
  106. %       an epsilon value for time interval calculation is created.
  107.  
  108. teps = intstep*1e-8;
  109.  
  110. %       Now check (or select) a final time and create the output
  111. %       time vector.  
  112.  
  113. if nargin < 3,
  114.     tfinal = max(ut);
  115. elseif tfinal < min(ut),
  116.     error('final time less than initial time')
  117.     return
  118.     end
  119. yt = [ut(1):intstep:tfinal].';
  120. nypts = length(yt);
  121.  
  122. %       Check if the final time is less than (to a tolerance of
  123. %       the integration step size) the maximum time in the input.
  124. %       If it is truncate the input.  Doing this now can save a
  125. %       lot of time in the interpolation.
  126.  
  127. if max(ut) - yt(nypts) > teps,
  128.     nupts = max(find(ut <= yt(nypts)+teps));
  129.     ut = ut(1:nupts);
  130.     udat = udat(1:uptr(nupts)+nru-1,:);
  131.     end
  132.  
  133. %       Now interpolate the input vector to the integration
  134. %       step size.   The only case where this is not necessary
  135. %       is when the stepsize is equal the spacing of a
  136. %       regular input vector.
  137.  
  138. interflg = 0;
  139. if nypts ~= nupts,
  140.     interflg = 1;
  141. elseif max(abs(yt - ut)) > teps,
  142.     interflg = 1;
  143.     end
  144.  
  145. if interflg == 1,
  146.     uint = zeros(nru*nypts,1);
  147.     disp('interpolating input vector (zero order hold)')
  148.     if nypts <= nupts,
  149.         for i = 1:nypts,
  150.             idx = max(find(ut <= yt(i)+teps));
  151.             uint((i-1)*nru+1:i*nru) = udat((idx-1)*nru+1:idx*nru);
  152.             end
  153.     elseif nupts == 1,
  154.         uint = reshape(udat*ones(1,nypts),nru*nypts,1);
  155.     else
  156.         for i = 1:nupts-1,
  157.             idx = find(yt >= ut(i) & yt < ut(i+1));
  158.             uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
  159.               reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
  160.               nru*length(idx),1);
  161.             end
  162.         idx = find(yt >= ut(nupts));
  163.         uint((min(idx)-1)*nru+1:max(idx)*nru) = ...
  164.           reshape(udat((i-1)*nru+1:i*nru)*ones(1,length(idx)), ...
  165.           nru*length(idx),1);
  166.         end
  167.     udat = uint;
  168.     end
  169.  
  170. %       reshape the input so that the matrix multiplications
  171. %       line up.  At this point yt should correspond to ut.
  172.  
  173. udat = reshape(udat,nru,nypts);
  174.     
  175. %       the input vector has now been expanded to a matrix: u
  176. %       constant systems are treated by a multiplication call
  177. %       For a,b,c,d quadruples a digitization is done.
  178.  
  179. y = zeros(nypts*ny+1,2);
  180. if type == 'syst',
  181.     x = [x0 zeros(length(x0),nypts-1)];
  182.     end
  183.  
  184. %       First, the single time point case is dealt with.  Here
  185. %       only the d term of a system and x0 affect the output.
  186.  
  187. if nypts == 1,
  188.     if type == 'syst',
  189.         y(1:ny,1) = c*x + d*udat;
  190.     else
  191.         y(1:ny,1) = sys*udat;
  192.         end
  193.  
  194. %       Now the multiple timepoint case is treated.  Here 
  195. %       a check is made to determine if the system has dynamics.
  196. %       The no dynamics case can be handled by a simple multiplication.
  197.  
  198. elseif type == 'syst',
  199.     ab = [a*intstep b*intstep;zeros(nu,nx+nu)];
  200.     eab = expm(ab);
  201.     ad = eab(1:nx,1:nx);
  202.     bd = eab(1:nx,nx+1:nx+nu);
  203.     for i=2:nypts,
  204.         x(:,i) = ad*x(:,i-1) + bd*udat(:,i-1);
  205.         end
  206.     y(1:ny*nypts,1) = reshape(c*x + d*udat,ny*nypts,1);
  207. else
  208.     y(1:ny*nypts,1) = reshape(sys*udat,ny*nypts,1);
  209.     end
  210.  
  211. %       The output is reshaped into the usual varying format and 
  212. %       the time is packed with it.  In the case of only a single
  213. %       time point, a varying matrix is still created.  This is
  214. %       after all, assumed to be a time response.
  215.  
  216. y(1:nypts,2) = yt;
  217. y((ny*nypts)+1,2) = inf;
  218. y((ny*nypts)+1,1) = nypts;
  219. %
  220. % Copyright MUSYN INC 1991,  All Rights Reserved
  221.