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