|
3431
|
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 |
|
7016
|
6 ## under the terms of the GNU General Public License as published by |
|
|
7 ## the Free Software Foundation; either version 3 of the License, or (at |
|
|
8 ## your option) any later version. |
|
3431
|
9 ## |
|
7016
|
10 ## Octave is distributed in the hope that it will be useful, but |
|
|
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
|
|
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
|
|
13 ## General Public License for more details. |
|
3431
|
14 ## |
|
|
15 ## You should have received a copy of the GNU General Public License |
|
7016
|
16 ## along with Octave; see the file COPYING. If not, see |
|
|
17 ## <http://www.gnu.org/licenses/>. |
|
3431
|
18 |
|
|
19 ## -*- texinfo -*- |
|
3500
|
20 ## @deftypefn {Function File} {} d2c (@var{sys}, @var{tol}) |
|
|
21 ## @deftypefnx {Function File} {} d2c (@var{sys}, @var{opt}) |
|
5016
|
22 ## Convert a discrete (sub)system into a purely continuous one. |
|
|
23 ## The sampling time used is @code{sysgettsam(@var{sys})}. |
|
3431
|
24 ## |
|
|
25 ## @strong{Inputs} |
|
|
26 ## @table @var |
|
|
27 ## @item sys |
|
|
28 ## system data structure with discrete components |
|
|
29 ## @item tol |
|
|
30 ## Scalar value. |
|
5016
|
31 ## Tolerance for convergence of default @code{"log"} option (see below) |
|
3431
|
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 |
|
3502
|
39 ## @var{a}, @var{b}, to get a better fit to the original data. |
|
3431
|
40 ## |
|
3500
|
41 ## If called as @code{d2c (@var{sys}, @var{tol})}, with @var{tol} |
|
|
42 ## positive scalar, the @code{"log"} option is used. The default value |
|
|
43 ## for @var{tol} is @code{1e-8}. |
|
3431
|
44 ## @item "bi" |
|
|
45 ## Conversion is performed via bilinear transform |
|
3502
|
46 ## @math{z = (1 + s T / 2)/(1 - s T / 2)} where @math{T} is the |
|
3431
|
47 ## system sampling time (see @code{sysgettsam}). |
|
|
48 ## |
|
|
49 ## FIXME: bilinear option exits with an error if @var{sys} is not purely |
|
|
50 ## discrete |
|
|
51 ## @end table |
|
|
52 ## @end table |
|
5016
|
53 ## @strong{Output} |
|
|
54 ## @table @var |
|
|
55 ## @item csys |
|
|
56 ## continuous time system (same dimensions and signal names as in @var{sys}). |
|
|
57 ## @end table |
|
3431
|
58 ## @end deftypefn |
|
|
59 |
|
|
60 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
61 ## Created: August 23, 1994 |
|
|
62 ## Updated by John Ingram for system data structure August 1996 |
|
|
63 |
|
|
64 function csys = d2c (sys, opt) |
|
|
65 |
|
|
66 ## SYS_INTERNAL accesses members of system data structure |
|
|
67 |
|
|
68 if( (nargin != 1) & (nargin != 2) ) |
|
6046
|
69 print_usage (); |
|
4030
|
70 elseif (!isstruct(sys)) |
|
3431
|
71 error("sys must be in system data structure"); |
|
|
72 elseif(nargin == 1) |
|
|
73 opt = "log"; |
|
|
74 tol = 1e-12; |
|
5443
|
75 elseif(ischar(opt)) # all remaining cases are for nargin == 2 |
|
3431
|
76 tol = 1e-12; |
|
|
77 if( !(strcmp(opt,"log") | strcmp(opt,"bi") ) ) |
|
|
78 error(["d2c: invalid opt passed=",opt]); |
|
|
79 endif |
|
|
80 elseif(!is_sample(opt)) |
|
7001
|
81 error("tol must be a positive scalar") |
|
3431
|
82 elseif(opt > 1e-2) |
|
|
83 warning(["d2c: ridiculous error tolerance passed=",num2str(opt); ... |
|
|
84 ", intended c2d call?"]) |
|
|
85 else |
|
|
86 tol = opt; |
|
|
87 opt = "log"; |
|
|
88 endif |
|
|
89 T = sysgettsam(sys); |
|
|
90 |
|
|
91 if(strcmp(opt,"bi")) |
|
|
92 ## bilinear transform |
|
|
93 ## convert with bilinear transform |
|
|
94 if (! is_digital(sys) ) |
|
|
95 error("d2c requires a discrete time system for input") |
|
|
96 endif |
|
|
97 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys); |
|
|
98 |
|
|
99 poles = eig(a); |
|
|
100 if( find(abs(poles-1) < 200*(n+nz)*eps) ) |
|
|
101 warning("d2c: some poles very close to one. May get bad results."); |
|
|
102 endif |
|
|
103 |
|
|
104 I = eye(size(a)); |
|
|
105 tk = 2/sqrt(T); |
|
|
106 A = (2/T)*(a-I)/(a+I); |
|
|
107 iab = (I+a)\b; |
|
|
108 B = tk*iab; |
|
|
109 C = tk*(c/(I+a)); |
|
|
110 D = d- (c*iab); |
|
|
111 stnamec = strappend(stname,"_c"); |
|
4771
|
112 csys = ss(A,B,C,D,0,rows(A),0,stnamec,inname,outname); |
|
3431
|
113 elseif(strcmp(opt,"log")) |
|
|
114 sys = sysupdate(sys,"ss"); |
|
|
115 [n,nz,m,p] = sysdimensions(sys); |
|
|
116 |
|
|
117 if(nz == 0) |
|
|
118 warning("d2c: all states continuous; setting outputs to agree"); |
|
|
119 csys = syssetsignals(sys,"yd",zeros(1,1:p)); |
|
|
120 return; |
|
|
121 elseif(n != 0) |
|
|
122 warning(["d2c: n=",num2str(n),">0; performing c2d first"]); |
|
|
123 sys = c2d(sys,T); |
|
|
124 endif |
|
|
125 [a,b] = sys2ss(sys); |
|
|
126 |
|
|
127 [ma,na] = size(a); |
|
|
128 [mb,nb] = size(b); |
|
|
129 |
|
|
130 if(isempty(b) ) |
|
|
131 warning("d2c: empty b matrix"); |
|
|
132 Amat = a; |
|
|
133 else |
|
|
134 Amat = [a, b; zeros(nb,na), eye(nb)]; |
|
|
135 endif |
|
|
136 |
|
|
137 poles = eig(a); |
|
|
138 if( find(abs(poles) < 200*(n+nz)*eps) ) |
|
|
139 warning("d2c: some poles very close to zero. logm not performed"); |
|
|
140 Mtop = zeros(ma, na+nb); |
|
|
141 elseif( find(abs(poles-1) < 200*(n+nz)*eps) ) |
|
|
142 warning("d2c: some poles very close to one. May get bad results."); |
|
|
143 logmat = real(logm(Amat)/T); |
|
|
144 Mtop = logmat(1:na,:); |
|
|
145 else |
|
|
146 logmat = real(logm(Amat)/T); |
|
|
147 Mtop = logmat(1:na,:); |
|
|
148 endif |
|
|
149 |
|
|
150 ## perform simplistic, stupid optimization approach. |
|
|
151 ## should re-write with a Davidson-Fletcher CG approach |
|
|
152 mxthresh = norm(Mtop); |
|
|
153 if(mxthresh == 0) |
|
|
154 mxthresh = 1; |
|
|
155 endif |
|
|
156 eps1 = mxthresh; #gradient descent step size |
|
|
157 cnt = max(20,(n*nz)*4); #max number of iterations |
|
|
158 newgrad=1; #signal for new gradient |
|
|
159 while( (eps1/mxthresh > tol) & cnt) |
|
|
160 cnt = cnt-1; |
|
|
161 ## calculate the gradient of error with respect to Amat... |
|
|
162 geps = norm(Mtop)*1e-8; |
|
|
163 if(geps == 0) |
|
|
164 geps = 1e-8; |
|
|
165 endif |
|
|
166 DMtop = Mtop; |
|
|
167 if(isempty(b)) |
|
|
168 Mall = Mtop; |
|
|
169 DMall = DMtop; |
|
|
170 else |
|
|
171 Mall = [Mtop; zeros(nb,na+nb)]; |
|
|
172 DMall = [DMtop; zeros(nb,na+nb) ]; |
|
|
173 endif |
|
|
174 |
|
|
175 if(newgrad) |
|
|
176 GrMall = zeros(size(Mall)); |
|
|
177 for ii=1:rows(Mtop) |
|
|
178 for jj=1:columns(Mtop) |
|
|
179 DMall(ii,jj) = Mall(ii,jj) + geps; |
|
|
180 GrMall(ii,jj) = norm (Amat - expm (DMall*T), "fro") ... |
|
|
181 - norm (Amat - expm (Mall*T), "fro"); |
|
|
182 DMall(ii,jj) = Mall(ii,jj); |
|
|
183 endfor |
|
|
184 endfor |
|
|
185 GrMall = GrMall/norm(GrMall,1); |
|
|
186 newgrad = 0; |
|
|
187 endif |
|
|
188 |
|
|
189 ## got a gradient, now try to use it |
|
|
190 DMall = Mall-eps1*GrMall; |
|
|
191 |
|
|
192 FMall = expm(Mall*T); |
|
|
193 FDMall = expm(DMall*T); |
|
|
194 FmallErr = norm(Amat - FMall); |
|
|
195 FdmallErr = norm(Amat - FDMall); |
|
|
196 if( FdmallErr < FmallErr) |
|
|
197 Mtop = DMall(1:na,:); |
|
|
198 eps1 = min(eps1*2,1e12); |
|
|
199 newgrad = 1; |
|
|
200 else |
|
|
201 eps1 = eps1/2; |
|
|
202 endif |
|
|
203 |
|
|
204 if(FmallErr == 0) |
|
|
205 eps1 = 0; |
|
|
206 endif |
|
|
207 |
|
|
208 endwhile |
|
|
209 |
|
|
210 [aa,bb,cc,dd,tsam,nn,nz,stnam,innam,outnam,yd] = sys2ss(sys); |
|
|
211 aa = Mall(1:na,1:na); |
|
|
212 if(!isempty(b)) |
|
|
213 bb = Mall(1:na,(na+1):(na+nb)); |
|
|
214 endif |
|
4771
|
215 csys = ss(aa,bb,cc,dd,0,na,0,stnam,innam,outnam); |
|
3431
|
216 |
|
|
217 ## update names |
|
|
218 nn = sysdimensions(sys); |
|
|
219 for ii = (nn+1):na |
|
|
220 strval = sprintf("%s_c",sysgetsignals(csys,"st",ii,1)); |
|
|
221 csys = syssetsignals(csys,"st",strval,ii); |
|
|
222 endfor |
|
|
223 endif |
|
|
224 |
|
|
225 endfunction |
