Mercurial > hg > octave-nkf

comparison scripts/control/d2c.m @ 3381:69b167451491

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1999-12-15 20:48:10 by jwe]
author jwe
date Wed, 15 Dec 1999 20:48:45 +0000
parents 8dd4718801fd
children ebf69eb3d07d
comparison
equal deleted inserted replaced
3380:f5edd74bdc6c 3381:69b167451491
1 # Copyright (C) 1996,1998 Auburn University. All Rights Reserved 1 ## Copyright (C) 1996,1998 Auburn University. All Rights Reserved
2 # 2 ##
3 # This file is part of Octave. 3 ## This file is part of Octave.
4 # 4 ##
5 # Octave is free software; you can redistribute it and/or modify it 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 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 7 ## Free Software Foundation; either version 2, or (at your option) any
8 # later version. 8 ## later version.
9 # 9 ##
10 # Octave is distributed in the hope that it will be useful, but WITHOUT 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 11 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 # FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 12 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 # for more details. 13 ## for more details.
14 # 14 ##
15 # You should have received a copy of the GNU General Public License 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 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. 17 ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
18 18
19 ## -*- texinfo -*- 19 ## -*- texinfo -*-
20 ## @deftypefn {Function File } {@var{csys} =} d2c (@var{sys}@{,@var{tol}@}) 20 ## @deftypefn {Function File } {@var{csys} =} d2c (@var{sys}@{,@var{tol}@})
21 ## @deftypefnx {Function File } {@var{csys} =} d2c (@var{sys}, @var{opt}) 21 ## @deftypefnx {Function File } {@var{csys} =} d2c (@var{sys}, @var{opt})
22 ## Convert discrete (sub)system to a purely continuous system. Sampling 22 ## Convert discrete (sub)system to a purely continuous system. Sampling
55 ## @end deftypefn 55 ## @end deftypefn
56 ## 56 ##
57 57
58 58
59 function csys = d2c(sys,opt) 59 function csys = d2c(sys,opt)
60 # Written by R. Bruce Tenison August 23, 1994 60
61 # Updated by John Ingram for system data structure August 1996 61 ## Written by R. Bruce Tenison August 23, 1994
62 # SYS_INTERNAL accesses members of system data structure 62 ## Updated by John Ingram for system data structure August 1996
63 ## SYS_INTERNAL accesses members of system data structure
63 64
64 save_val = implicit_str_to_num_ok; # save for later 65 save_val = implicit_str_to_num_ok; # save for later
65 implicit_str_to_num_ok = 1; 66 implicit_str_to_num_ok = 1;
66 67
67 if( (nargin != 1) & (nargin != 2) ) 68 if( (nargin != 1) & (nargin != 2) )
86 opt = "log"; 87 opt = "log";
87 endif 88 endif
88 T = sysgettsam(sys); 89 T = sysgettsam(sys);
89 90
90 if(strcmp(opt,"bi")) 91 if(strcmp(opt,"bi"))
91 # bilinear transform 92 ## bilinear transform
92 # convert with bilinear transform 93 ## convert with bilinear transform
93 if (! is_digital(sys) ) 94 if (! is_digital(sys) )
94 error("d2c requires a discrete time system for input") 95 error("d2c requires a discrete time system for input")
95 endif 96 endif
96 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys); 97 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys);
97 98
144 else 145 else
145 logmat = real(logm(Amat)/T); 146 logmat = real(logm(Amat)/T);
146 Mtop = logmat(1:na,:); 147 Mtop = logmat(1:na,:);
147 endif 148 endif
148 149
149 # perform simplistic, stupid optimization approach. 150 ## perform simplistic, stupid optimization approach.
150 # should re-write with a Davidson-Fletcher CG approach 151 ## should re-write with a Davidson-Fletcher CG approach
151 mxthresh = norm(Mtop); 152 mxthresh = norm(Mtop);
152 if(mxthresh == 0) 153 if(mxthresh == 0)
153 mxthresh = 1; 154 mxthresh = 1;
154 endif 155 endif
155 eps1 = mxthresh; #gradient descent step size 156 eps1 = mxthresh; #gradient descent step size
156 cnt = max(20,(n*nz)*4); #max number of iterations 157 cnt = max(20,(n*nz)*4); #max number of iterations
157 newgrad=1; #signal for new gradient 158 newgrad=1; #signal for new gradient
158 while( (eps1/mxthresh > tol) & cnt) 159 while( (eps1/mxthresh > tol) & cnt)
159 cnt = cnt-1; 160 cnt = cnt-1;
160 # calculate the gradient of error with respect to Amat... 161 ## calculate the gradient of error with respect to Amat...
161 geps = norm(Mtop)*1e-8; 162 geps = norm(Mtop)*1e-8;
162 if(geps == 0) 163 if(geps == 0)
163 geps = 1e-8; 164 geps = 1e-8;
164 endif 165 endif
165 DMtop = Mtop; 166 DMtop = Mtop;
183 endfor 184 endfor
184 GrMall = GrMall/norm(GrMall,1); 185 GrMall = GrMall/norm(GrMall,1);
185 newgrad = 0; 186 newgrad = 0;
186 endif 187 endif
187 188
188 #got a gradient, now try to use it 189 ## got a gradient, now try to use it
189 DMall = Mall-eps1*GrMall; 190 DMall = Mall-eps1*GrMall;
190 191
191 FMall = expm(Mall*T); 192 FMall = expm(Mall*T);
192 FDMall = expm(DMall*T); 193 FDMall = expm(DMall*T);
193 FmallErr = norm(Amat - FMall); 194 FmallErr = norm(Amat - FMall);
211 if(!isempty(b)) 212 if(!isempty(b))
212 bb = Mall(1:na,(na+1):(na+nb)); 213 bb = Mall(1:na,(na+1):(na+nb));
213 endif 214 endif
214 csys = ss2sys(aa,bb,cc,dd,0,na,0,stnam,innam,outnam); 215 csys = ss2sys(aa,bb,cc,dd,0,na,0,stnam,innam,outnam);
215 216
216 # update names 217 ## update names
217 nn = sysdimensions(sys); 218 nn = sysdimensions(sys);
218 for ii = (nn+1):na 219 for ii = (nn+1):na
219 strval = sprintf("%s_c",sysgetsignals(csys,"st",ii,1)); 220 strval = sprintf("%s_c",sysgetsignals(csys,"st",ii,1));
220 csys = syssetsignals(csys,"st",strval,ii); 221 csys = syssetsignals(csys,"st",strval,ii);
221 endfor 222 endfor