|
3279
|
1 # Copyright (C) 1996,1998 Auburn University. All Rights Reserved |
|
3213
|
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 |
|
3284
|
17 # Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. |
|
3346
|
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{sysgettsam(@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{sysgettsam}). |
|
|
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 |
|
3213
|
58 |
|
|
59 function csys = d2c(sys,opt) |
|
|
60 # Written by R. Bruce Tenison August 23, 1994 |
|
|
61 # Updated by John Ingram for system data structure August 1996 |
|
|
62 # SYS_INTERNAL accesses members of system data structure |
|
|
63 |
|
|
64 save_val = implicit_str_to_num_ok; # save for later |
|
|
65 implicit_str_to_num_ok = 1; |
|
|
66 |
|
|
67 if( (nargin != 1) & (nargin != 2) ) |
|
|
68 usage("csys = d2c(sys[,tol]), csys = d2c(sys,opt)"); |
|
|
69 elseif (!is_struct(sys)) |
|
|
70 error("sys must be in system data structure"); |
|
|
71 elseif(nargin == 1) |
|
|
72 opt = "log"; |
|
|
73 tol = 1e-12; |
|
|
74 elseif(isstr(opt)) # all remaining cases are for nargin == 2 |
|
|
75 tol = 1e-12; |
|
|
76 if( !(strcmp(opt,"log") | strcmp(opt,"bi") ) ) |
|
|
77 error(["d2c: illegal opt passed=",opt]); |
|
|
78 endif |
|
|
79 elseif(!is_sample(opt)) |
|
|
80 error("tol must be a postive scalar") |
|
|
81 elseif(opt > 1e-2) |
|
|
82 warning(["d2c: ridiculous error tolerance passed=",num2str(opt); ... |
|
|
83 ", intended c2d call?"]) |
|
|
84 else |
|
|
85 tol = opt; |
|
|
86 opt = "log"; |
|
|
87 endif |
|
3228
|
88 T = sysgettsam(sys); |
|
3213
|
89 |
|
3228
|
90 if(strcmp(opt,"bi")) |
|
3213
|
91 # bilinear transform |
|
|
92 # convert with bilinear transform |
|
|
93 if (! is_digital(sys) ) |
|
|
94 error("d2c requires a discrete time system for input") |
|
|
95 endif |
|
|
96 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys); |
|
|
97 |
|
|
98 poles = eig(a); |
|
|
99 if( find(abs(poles-1) < 200*(n+nz)*eps) ) |
|
|
100 warning("d2c: some poles very close to one. May get bad results."); |
|
|
101 endif |
|
|
102 |
|
|
103 I = eye(size(a)); |
|
|
104 tk = 2/sqrt(T); |
|
|
105 A = (2/T)*(a-I)/(a+I); |
|
|
106 iab = (I+a)\b; |
|
|
107 B = tk*iab; |
|
|
108 C = tk*(c/(I+a)); |
|
|
109 D = d- (c*iab); |
|
3228
|
110 stnamec = strappend(stname,"_c"); |
|
3213
|
111 csys = ss2sys(A,B,C,D,0,rows(A),0,stnamec,inname,outname); |
|
3228
|
112 elseif(strcmp(opt,"log")) |
|
3213
|
113 sys = sysupdate(sys,"ss"); |
|
|
114 [n,nz,m,p] = sysdimensions(sys); |
|
|
115 |
|
|
116 if(nz == 0) |
|
|
117 warning("d2c: all states continuous; setting outputs to agree"); |
|
3228
|
118 csys = syssetsignals(sys,"yd",zeros(1,1:p)); |
|
3213
|
119 return; |
|
|
120 elseif(n != 0) |
|
|
121 warning(["d2c: n=",num2str(n),">0; performing c2d first"]); |
|
|
122 sys = c2d(sys,T); |
|
|
123 endif |
|
|
124 [a,b] = sys2ss(sys); |
|
|
125 |
|
|
126 [ma,na] = size(a); |
|
|
127 [mb,nb] = size(b); |
|
|
128 |
|
|
129 if(isempty(b) ) |
|
|
130 warning("d2c: empty b matrix"); |
|
|
131 Amat = a; |
|
|
132 else |
|
3238
|
133 Amat = [a, b; zeros(nb,na), eye(nb)]; |
|
3213
|
134 endif |
|
|
135 |
|
|
136 poles = eig(a); |
|
|
137 if( find(abs(poles) < 200*(n+nz)*eps) ) |
|
|
138 warning("d2c: some poles very close to zero. logm not performed"); |
|
|
139 Mtop = zeros(ma, na+nb); |
|
|
140 elseif( find(abs(poles-1) < 200*(n+nz)*eps) ) |
|
|
141 warning("d2c: some poles very close to one. May get bad results."); |
|
|
142 logmat = real(logm(Amat)/T); |
|
|
143 Mtop = logmat(1:na,:); |
|
|
144 else |
|
|
145 logmat = real(logm(Amat)/T); |
|
|
146 Mtop = logmat(1:na,:); |
|
|
147 endif |
|
|
148 |
|
|
149 # perform simplistic, stupid optimization approach. |
|
|
150 # should re-write with a Davidson-Fletcher CG approach |
|
|
151 mxthresh = norm(Mtop); |
|
|
152 if(mxthresh == 0) |
|
|
153 mxthresh = 1; |
|
|
154 endif |
|
|
155 eps1 = mxthresh; #gradient descent step size |
|
|
156 cnt = max(20,(n*nz)*4); #max number of iterations |
|
|
157 newgrad=1; #signal for new gradient |
|
|
158 while( (eps1/mxthresh > tol) & cnt) |
|
|
159 cnt = cnt-1; |
|
|
160 # calculate the gradient of error with respect to Amat... |
|
|
161 geps = norm(Mtop)*1e-8; |
|
|
162 if(geps == 0) |
|
|
163 geps = 1e-8; |
|
|
164 endif |
|
|
165 DMtop = Mtop; |
|
|
166 if(isempty(b)) |
|
|
167 Mall = Mtop; |
|
|
168 DMall = DMtop; |
|
|
169 else |
|
3238
|
170 Mall = [Mtop; zeros(nb,na+nb)]; |
|
|
171 DMall = [DMtop; zeros(nb,na+nb) ]; |
|
3213
|
172 endif |
|
|
173 |
|
|
174 if(newgrad) |
|
|
175 GrMall = zeros(size(Mall)); |
|
|
176 for ii=1:rows(Mtop) |
|
|
177 for jj=1:columns(Mtop) |
|
|
178 DMall(ii,jj) = Mall(ii,jj) + geps; |
|
|
179 GrMall(ii,jj) = norm(Amat - expm(DMall*T),'fro') ... |
|
|
180 - norm(Amat-expm(Mall*T),'fro'); |
|
|
181 DMall(ii,jj) = Mall(ii,jj); |
|
|
182 endfor |
|
|
183 endfor |
|
|
184 GrMall = GrMall/norm(GrMall,1); |
|
|
185 newgrad = 0; |
|
|
186 endif |
|
|
187 |
|
|
188 #got a gradient, now try to use it |
|
|
189 DMall = Mall-eps1*GrMall; |
|
|
190 |
|
|
191 FMall = expm(Mall*T); |
|
|
192 FDMall = expm(DMall*T); |
|
|
193 FmallErr = norm(Amat - FMall); |
|
|
194 FdmallErr = norm(Amat - FDMall); |
|
|
195 if( FdmallErr < FmallErr) |
|
|
196 Mtop = DMall(1:na,:); |
|
|
197 eps1 = min(eps1*2,1e12); |
|
|
198 newgrad = 1; |
|
|
199 else |
|
|
200 eps1 = eps1/2; |
|
|
201 endif |
|
|
202 |
|
|
203 if(FmallErr == 0) |
|
|
204 eps1 = 0; |
|
|
205 endif |
|
|
206 |
|
|
207 endwhile |
|
|
208 |
|
3228
|
209 [aa,bb,cc,dd,tsam,nn,nz,stnam,innam,outnam,yd] = sys2ss(sys); |
|
|
210 aa = Mall(1:na,1:na); |
|
3213
|
211 if(!isempty(b)) |
|
3228
|
212 bb = Mall(1:na,(na+1):(na+nb)); |
|
3213
|
213 endif |
|
3228
|
214 csys = ss2sys(aa,bb,cc,dd,0,na,0,stnam,innam,outnam); |
|
3213
|
215 |
|
3228
|
216 # update names |
|
|
217 nn = sysdimensions(sys); |
|
|
218 for ii = (nn+1):na |
|
|
219 strval = sprintf("%s_c",sysgetsignals(csys,"st",ii,1)); |
|
|
220 csys = syssetsignals(csys,"st",strval,ii); |
|
3213
|
221 endfor |
|
|
222 endif |
|
|
223 |
|
|
224 implicit_str_to_num_ok = save_val; # restore value |
|
|
225 endfunction |
