home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts.fat / control / d2c.m < prev    next >
Text File  |  1999-12-24  |  7KB  |  223 lines
  1. ## Copyright (C) 1996,1998 Auburn University.  All Rights Reserved
  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, 59 Temple Place, Suite 330, Boston, MA 02111 USA. 
  18.  
  19. ## -*- texinfo -*-
  20. ## @deftypefn {Function File } {@var{csys} =} d2c (@var{sys}@{,@var{tol}@})
  21. ## @deftypefnx {Function File } {@var{csys} =} d2c (@var{sys}, @var{opt})
  22. ## Convert discrete (sub)system to a purely continuous system.  Sampling
  23. ## time used is @code{sysgetts(@var{sys})}
  24. ## 
  25. ## @strong{Inputs}
  26. ## @table @var
  27. ## @item   sys
  28. ##  system data structure with discrete components
  29. ## @item   tol
  30. ## Scalar value.
  31. ##  tolerance for convergence of default @code{"log"} option (see below)
  32. ## @item   opt
  33. ##  conversion option.  Choose from:
  34. ## @table @code
  35. ## @item         "log"
  36. ##  (default) Conversion is performed via a matrix logarithm.
  37. ## Due to some problems with this computation, it is
  38. ## followed by a steepest descent algorithm to identify continuous time 
  39. ## @var{A}, @var{B}, to get a better fit to the original data.  
  40. ## 
  41. ## If called as @code{d2c}(@var{sys},@var{tol}), @var{tol=}positive scalar, 
  42. ##     the @code{"log"} option is used.  The default value for @var{tol} is 
  43. ##     @code{1e-8}.
  44. ## @item        "bi"
  45. ##  Conversion is performed via bilinear transform 
  46. ## @math{z = (1 + s T / 2)/(1 - s T / 2)} where @var{T} is the 
  47. ## system sampling time (see @code{sysgetts}).
  48. ## 
  49. ## FIXME: bilinear option exits with an error if @var{sys} is not purely discrete
  50. ## 
  51. ## @end table
  52. ## @end table
  53. ## @strong{Outputs} @var{csys} continuous time system (same dimensions and
  54. ## signal names as in @var{sys}).
  55. ## @end deftypefn
  56. ## 
  57.  
  58.  
  59. function csys = d2c (sys, opt)
  60.  
  61.   ## Written by R. Bruce Tenison August 23, 1994
  62.   ## Updated by John Ingram for system data structure  August 1996
  63.   ## SYS_INTERNAL accesses members of system data structure
  64.  
  65.   if( (nargin != 1) & (nargin != 2) )
  66.     usage("csys = d2c(sys[,tol]), csys = d2c(sys,opt)");
  67.   elseif (!is_struct(sys))
  68.     error("sys must be in system data structure");
  69.   elseif(nargin == 1)
  70.     opt = "log";
  71.     tol = 1e-12;
  72.   elseif(isstr(opt))   # all remaining cases are for nargin == 2
  73.     tol = 1e-12;
  74.     if( !(strcmp(opt,"log") | strcmp(opt,"bi") ) )
  75.       error(["d2c: illegal opt passed=",opt]);
  76.     endif
  77.   elseif(!is_sampl(opt))
  78.     error("tol must be a postive scalar")
  79.   elseif(opt > 1e-2)
  80.     warning(["d2c: ridiculous error tolerance passed=",num2str(opt); ...
  81.     ", intended c2d call?"])
  82.   else
  83.     tol = opt;
  84.     opt = "log";
  85.   endif
  86.   T = sysgetts(sys);
  87.  
  88.   if(strcmp(opt,"bi"))
  89.     ## bilinear transform
  90.     ## convert with bilinear transform
  91.     if (! is_digit(sys) )
  92.        error("d2c requires a discrete time system for input")
  93.     endif
  94.     [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys);
  95.  
  96.     poles = eig(a);
  97.     if( find(abs(poles-1) < 200*(n+nz)*eps) )
  98.       warning("d2c: some poles very close to one.  May get bad results.");
  99.     endif
  100.  
  101.     I = eye(size(a));
  102.     tk = 2/sqrt(T);
  103.     A = (2/T)*(a-I)/(a+I);
  104.     iab = (I+a)\b;
  105.     B = tk*iab;
  106.     C = tk*(c/(I+a));
  107.     D = d- (c*iab);
  108.     stnamec = strappen(stname,"_c");
  109.     csys = ss2sys(A,B,C,D,0,rows(A),0,stnamec,inname,outname);
  110.   elseif(strcmp(opt,"log"))
  111.     sys = sysupdat(sys,"ss");
  112.     [n,nz,m,p] = sysdimen(sys);
  113.   
  114.     if(nz == 0)
  115.       warning("d2c: all states continuous; setting outputs to agree");
  116.       csys = syssetsg(sys,"yd",zeros(1,1:p));
  117.       return;
  118.     elseif(n != 0)
  119.       warning(["d2c: n=",num2str(n),">0; performing c2d first"]);
  120.       sys = c2d(sys,T);
  121.     endif
  122.     [a,b] = sys2ss(sys);
  123.   
  124.     [ma,na] = size(a);
  125.     [mb,nb] = size(b);
  126.   
  127.     if(isempty(b) )
  128.       warning("d2c: empty b matrix");
  129.       Amat = a;
  130.     else
  131.       Amat = [a, b; zeros(nb,na), eye(nb)];
  132.     endif
  133.   
  134.     poles = eig(a);
  135.     if( find(abs(poles) < 200*(n+nz)*eps) )
  136.       warning("d2c: some poles very close to zero.  logm not performed");
  137.       Mtop = zeros(ma, na+nb);
  138.     elseif( find(abs(poles-1) < 200*(n+nz)*eps) )
  139.       warning("d2c: some poles very close to one.  May get bad results.");
  140.       logmat = real(logm(Amat)/T);
  141.       Mtop = logmat(1:na,:);
  142.     else
  143.       logmat = real(logm(Amat)/T);
  144.       Mtop = logmat(1:na,:);
  145.     endif
  146.   
  147.     ## perform simplistic, stupid optimization approach.
  148.     ## should re-write with a Davidson-Fletcher CG approach
  149.     mxthresh = norm(Mtop);
  150.     if(mxthresh == 0)
  151.       mxthresh = 1;
  152.     endif
  153.     eps1 = mxthresh;    #gradient descent step size
  154.     cnt = max(20,(n*nz)*4);    #max number of iterations
  155.     newgrad=1;    #signal for new gradient
  156.     while( (eps1/mxthresh > tol) & cnt)
  157.       cnt = cnt-1;
  158.       ## calculate the gradient of error with respect to Amat...
  159.       geps = norm(Mtop)*1e-8;
  160.       if(geps == 0)
  161.         geps = 1e-8;
  162.       endif
  163.       DMtop = Mtop;
  164.       if(isempty(b))
  165.         Mall = Mtop;
  166.         DMall = DMtop;
  167.       else
  168.         Mall = [Mtop; zeros(nb,na+nb)];
  169.         DMall = [DMtop; zeros(nb,na+nb) ];
  170.       endif
  171.   
  172.       if(newgrad)
  173.         GrMall = zeros(size(Mall));
  174.         for ii=1:rows(Mtop)
  175.           for jj=1:columns(Mtop)
  176.         DMall(ii,jj) = Mall(ii,jj) + geps;
  177.             GrMall(ii,jj) = norm(Amat - expm(DMall*T),'fro') ...
  178.           - norm(Amat-expm(Mall*T),'fro');
  179.             DMall(ii,jj) = Mall(ii,jj);
  180.           endfor
  181.         endfor
  182.         GrMall = GrMall/norm(GrMall,1);
  183.         newgrad = 0;
  184.       endif
  185.   
  186.       ## got a gradient, now try to use it
  187.       DMall = Mall-eps1*GrMall;
  188.   
  189.       FMall = expm(Mall*T);
  190.       FDMall = expm(DMall*T);
  191.       FmallErr = norm(Amat - FMall);
  192.       FdmallErr = norm(Amat - FDMall);
  193.       if( FdmallErr < FmallErr)
  194.         Mtop = DMall(1:na,:);
  195.         eps1 = min(eps1*2,1e12);
  196.         newgrad = 1;
  197.       else
  198.         eps1 = eps1/2;
  199.       endif
  200.   
  201.       if(FmallErr == 0)
  202.         eps1 = 0;
  203.       endif
  204.       
  205.     endwhile
  206.   
  207.     [aa,bb,cc,dd,tsam,nn,nz,stnam,innam,outnam,yd] = sys2ss(sys);
  208.     aa = Mall(1:na,1:na);
  209.     if(!isempty(b))
  210.       bb = Mall(1:na,(na+1):(na+nb));
  211.     endif
  212.     csys = ss2sys(aa,bb,cc,dd,0,na,0,stnam,innam,outnam);
  213.     
  214.     ## update names
  215.     nn = sysdimen(sys);
  216.     for ii = (nn+1):na
  217.       strval = sprintf("%s_c",sysgetsg(csys,"st",ii,1));
  218.       csys = syssetsg(csys,"st",strval,ii);
  219.     endfor
  220.   endif
  221.  
  222. endfunction
  223.