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