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