home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / control / c2d.m < prev    next >
Text File  |  1999-03-05  |  5KB  |  153 lines
  1. # Copyright (C) 1993, 1994, 1995 John W. Eaton
  2. # 
  3. # This file is part of Octave.
  4. # 
  5. # Octave is free software; you can redistribute it and/or modify it
  6. # under the terms of the GNU General Public License as published by the
  7. # Free Software Foundation; either version 2, or (at your option) any
  8. # later version.
  9. # 
  10. # Octave is distributed in the hope that it will be useful, but WITHOUT
  11. # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12. # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13. # for more details.
  14. # 
  15. # You should have received a copy of the GNU General Public License
  16. # along with Octave; see the file COPYING.  If not, write to the Free
  17. # Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  
  19. function dsys = c2d (sys, opt, T)
  20.  
  21. # Usage: dsys = c2d (sys[, T])
  22. # Usage: dsys = c2d (sys[, opt[, T]])
  23. # inputs:
  24. #   sys: system data structure (may be mixed discrete/continiuous time)
  25. #   optional arguments:
  26. #     opt: conversion option: 
  27. #          "ex" - use the matrix exponential (default)
  28. #          "bi" - use the bilinear transformation
  29. #                   2(z-1)
  30. #               s = -----
  31. #                   T(z+1)
  32. #               FIXME: This option exits with an error if sys is not purely 
  33. #               continuous. (The ex can handle mixed systems.)
  34. #
  35. #     T: sampling time; required if sys is purely continuous.
  36. # outputs: 
  37. #   dsys: discrete time equivalent via zero-order hold, sample each T sec.
  38. #
  39. # converts the system described by:
  40. #   .
  41. #   x = Ac x + Bc u
  42. #
  43. # into a discrete time equivalent model via the matrix exponential
  44. #
  45. #   x[n+1] = Ad x[n] + Bd u[n]
  46. #
  47. # Note: This function adds _d to the names of the new discrete states.   
  48.  
  49. # Written by R.B. Tenison (btenison@eng.auburn.edu)
  50. # October 1993
  51. # Updated by John Ingram for system data structure August 1996
  52.  
  53.   save_val = implicit_str_to_num_ok;    # save for later
  54.   implicit_str_to_num_ok = 1;
  55.  
  56. # parse input arguments
  57.   if(nargin < 1 | nargin > 3)
  58.     usage("dsys=c2d(sys[,T])");
  59.   elseif (!is_struct(sys))
  60.     error("sys must be a system data structure");
  61.   elseif (nargin == 1)
  62.     opt = "ex";
  63.   elseif (nargin == 2 & !isstr(opt) )
  64.     T = opt;
  65.     opt = "ex";
  66.   endif
  67.  
  68.   # check if sampling period T was passed.
  69.   Ts = sysgettsam(sys);
  70.   if(!exist("T"))
  71.     T = Ts;
  72.     if(T == 0)
  73.       error("sys is purely continuous; no sampling period T provided");
  74.     endif
  75.   elseif (T != Ts & Ts > 0)
  76.     warning(["c2d: T=",num2str(T),", system tsam==",num2str(Ts), ...
  77.       ": using T=", num2str(min(T,Ts))]);
  78.     T = min(T,Ts);
  79.   endif
  80.  
  81.   if (!is_sample(T))
  82.     error("sampling period T must be a postive, real scalar");
  83.   elseif( ! (strcmp(opt,"ex") | strcmp(opt,"bi") ) )
  84.     error(["illegal option passed: ",opt])
  85.   endif
  86.  
  87.   sys = sysupdate(sys,"ss");
  88.   [n,nz,m,p] = sysdimensions(sys);
  89.  
  90.   if (n == 0)
  91.     warning("c2d: sys has no continuous states; setting outputs to discrete");
  92.     dsys = syssetsignals(sys,"yd",ones(1:p));
  93.   elseif(strcmp(opt,"ex"))
  94.     # construct new state-space (a,b,c,d) for continuous subsystem
  95.     [csys,Acd] = syscont(sys);       # extract continuous subsystem
  96.     [csys_a, csys_b, csys_c, csys_d] = sys2ss(csys);
  97.     [ sys_a,  sys_b,  sys_c,  sys_d] = sys2ss( sys);
  98.     if(isempty(Acd))                Bmat = sys_b;
  99.     elseif(isempty(csys_b))         Bmat = Acd;
  100.     else                            Bmat = [Acd, csys_b];     endif
  101.     
  102.     row_zer = columns(Bmat);
  103.     csysn = sysdimensions(csys);
  104.     col_zer = csysn + row_zer;
  105.  
  106.     [csysa,csysb,csysc,csysd] = sys2ss(csys);
  107.     if(isempty(Bmat) )
  108.       warning("c2d: no inputs to continuous subsystem.");
  109.       mat = csysa;
  110.     else
  111.       mat = [csysa, Bmat ; zeros( row_zer,col_zer) ];
  112.     endif
  113.  
  114.     matexp = expm(mat * T);
  115.   
  116.     Abar = matexp( 1:csysn , 1:(csysn + columns(Acd)) );  
  117.     Bbar = matexp( 1:csysn , (columns(Abar) + 1):columns(matexp) );
  118.  
  119.     newnz = rows(Abar);
  120.     outlist = ones(1,rows(csysc));
  121.     [stnames,innames,outnames] = sysgetsignals(csys);
  122.     dsys = ss2sys(Abar,Bbar,csysc,csysd,T,0,newnz,stnames,innames, ...
  123.     outnames,outlist);
  124.     # rename states
  125.     for ii=1:newnz
  126.       strval = sprintf("%s_d",sysgetsignals(dsys,"st",ii,1));
  127.       dsys = syssetsignals(dsys,"st",strval,ii);
  128.     endfor
  129.  
  130.   elseif(strcmp(opt,"bi"))
  131.     if(is_digital(sys))
  132.       error("c2d: system is already digital")
  133.     else
  134.       # convert with bilinear transform
  135.       [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys);
  136.       IT = (2/T)*eye(size(a));
  137.       A = (IT+a)/(IT-a);
  138.       iab = (IT-a)\b;
  139.       tk=2/sqrt(T);
  140.       B = tk*iab;
  141.       C = tk*(c/(IT-a));
  142.       D = d + (c*iab);
  143.       stnamed = strappend(stname,"_d");
  144.       dsys = ss2sys(A,B,C,D,T,0,rows(A),stnamed,inname,outname);
  145.     endif
  146.   else
  147.     error(["Bad option=",opt])
  148.   endif
  149.   
  150.   implicit_str_to_num_ok = save_val;    # restore value
  151.  
  152. endfunction
  153.