|
3431
|
1 ## Copyright (C) 1993, 1994, 1995 John W. Eaton |
|
|
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 -*- |
|
3502
|
20 ## @deftypefn {Function File} {} c2d (@var{sys}, @var{opt}, @var{t}) |
|
|
21 ## @deftypefnx {Function File} {} c2d (@var{sys}, @var{t}) |
|
3431
|
22 ## |
|
5016
|
23 ## Converts the system data structure describing: |
|
|
24 ## @iftex |
|
|
25 ## @tex |
|
|
26 ## $$ \dot x = A_cx + B_cu $$ |
|
|
27 ## @end tex |
|
|
28 ## @end iftex |
|
|
29 ## @ifinfo |
|
|
30 ## @example |
|
|
31 ## . |
|
|
32 ## x = Ac x + Bc u |
|
|
33 ## @end example |
|
|
34 ## @end ifinfo |
|
|
35 ## into a discrete time equivalent model: |
|
|
36 ## @iftex |
|
|
37 ## @tex |
|
|
38 ## $$ x_{n+1} = A_dx_n + B_du_n $$ |
|
|
39 ## @end tex |
|
|
40 ## @end iftex |
|
|
41 ## @ifinfo |
|
|
42 ## @example |
|
|
43 ## x[n+1] = Ad x[n] + Bd u[n] |
|
|
44 ## @end example |
|
|
45 ## @end ifinfo |
|
|
46 ## via the matrix exponential or bilinear transform. |
|
|
47 ## |
|
3431
|
48 ## @strong{Inputs} |
|
|
49 ## @table @var |
|
|
50 ## @item sys |
|
|
51 ## system data structure (may have both continuous time and discrete |
|
|
52 ## time subsystems) |
|
|
53 ## @item opt |
|
|
54 ## string argument; conversion option (optional argument; |
|
|
55 ## may be omitted as shown above) |
|
|
56 ## @table @code |
|
|
57 ## @item "ex" |
|
|
58 ## use the matrix exponential (default) |
|
|
59 ## @item "bi" |
|
|
60 ## use the bilinear transformation |
|
5016
|
61 ## @iftex |
|
|
62 ## @tex |
|
|
63 ## $$ s = { 2(z-1) \over T(z+1) } $$ |
|
|
64 ## @end tex |
|
|
65 ## @end iftex |
|
|
66 ## @ifinfo |
|
3431
|
67 ## @example |
|
|
68 ## 2(z-1) |
|
|
69 ## s = ----- |
|
|
70 ## T(z+1) |
|
|
71 ## @end example |
|
5016
|
72 ## @end ifinfo |
|
3431
|
73 ## FIXME: This option exits with an error if @var{sys} is not purely |
|
|
74 ## continuous. (The @code{ex} option can handle mixed systems.) |
|
3759
|
75 ## @item "matched" |
|
|
76 ## Use the matched pole/zero equivalent transformation (currently only |
|
5016
|
77 ## works for purely continuous @acronym{SISO} systems). |
|
|
78 ## @end table |
|
|
79 ## @item t |
|
|
80 ## sampling time; required if @var{sys} is purely continuous. |
|
|
81 ## |
|
5479
|
82 ## @strong{Note} that if the second argument is not a string, @code{c2d()} |
|
5016
|
83 ## assumes that the second argument is @var{t} and performs |
|
|
84 ## appropriate argument checks. |
|
3431
|
85 ## @end table |
|
|
86 ## |
|
5016
|
87 ## @strong{Output} |
|
|
88 ## @table @var |
|
|
89 ## @item dsys |
|
|
90 ## Discrete time equivalent via zero-order hold, sample each @var{t} sec. |
|
|
91 ## @end table |
|
3431
|
92 ## |
|
4946
|
93 ## This function adds the suffix @code{_d} |
|
3431
|
94 ## to the names of the new discrete states. |
|
|
95 ## @end deftypefn |
|
|
96 |
|
|
97 ## Author: R. Bruce Tenison <btenison@eng.auburn.edu> |
|
|
98 ## Created: October 1993 |
|
|
99 ## Updated by John Ingram for system data structure August 1996 |
|
|
100 |
|
|
101 function dsys = c2d (sys, opt, T) |
|
|
102 ## parse input arguments |
|
|
103 if(nargin < 1 | nargin > 3) |
|
6046
|
104 print_usage (); |
|
4030
|
105 elseif (!isstruct(sys)) |
|
3431
|
106 error("sys must be a system data structure"); |
|
|
107 elseif (nargin == 1) |
|
|
108 opt = "ex"; |
|
5443
|
109 elseif (nargin == 2 & !ischar(opt) ) |
|
3431
|
110 T = opt; |
|
|
111 opt = "ex"; |
|
|
112 endif |
|
|
113 |
|
5443
|
114 if (! ischar (opt)) |
|
3759
|
115 error ("expecting option as a string"); |
|
|
116 endif |
|
|
117 |
|
3431
|
118 ## check if sampling period T was passed. |
|
|
119 Ts = sysgettsam(sys); |
|
|
120 if(!exist("T")) |
|
|
121 T = Ts; |
|
|
122 if(T == 0) |
|
|
123 error("sys is purely continuous; no sampling period T provided"); |
|
|
124 endif |
|
|
125 elseif (T != Ts & Ts > 0) |
|
|
126 warning(["c2d: T=",num2str(T),", system tsam==",num2str(Ts), ... |
|
|
127 ": using T=", num2str(min(T,Ts))]); |
|
|
128 T = min(T,Ts); |
|
|
129 endif |
|
|
130 |
|
|
131 if (!is_sample(T)) |
|
7001
|
132 error("sampling period T must be a positive, real scalar"); |
|
3759
|
133 elseif (! (strcmp (opt, "ex") |
|
|
134 || strcmp (opt, "bi") |
|
|
135 || strcmp (opt, "matched"))) |
|
|
136 error ("invalid option passed: %s", opt); |
|
3431
|
137 endif |
|
|
138 |
|
|
139 sys = sysupdate(sys,"ss"); |
|
|
140 [n,nz,m,p] = sysdimensions(sys); |
|
|
141 if(n == 0) |
|
|
142 dsys = syssetsignals(sys,"yd",ones(1:p)); |
|
|
143 elseif(strcmp(opt,"ex")); |
|
|
144 [aa,bb,cc,dd] = sys2ss(sys); |
|
|
145 crng= 1:n; |
|
|
146 drng = n+(1:nz); |
|
|
147 |
|
|
148 ## partition state equations into continuous, imaginary subsystems |
|
|
149 Ac = aa(crng,crng); |
|
|
150 Bc = bb(crng,:); |
|
|
151 if(nz == 0) |
|
|
152 Acd = Adc = Add = Bd = 0; |
|
|
153 else |
|
|
154 Acd = aa(crng,drng); |
|
|
155 Adc = aa(drng,crng); |
|
|
156 Add = aa(drng,drng); |
|
|
157 Bd = bb(drng,:); |
|
|
158 Bc = [Bc, Acd]; ## append discrete states as inputs to cont system |
|
|
159 endif |
|
|
160 |
|
|
161 ## convert state equations |
|
|
162 mat = [Ac, Bc; zeros(m+nz,n+nz+m)]; |
|
|
163 matexp = expm(mat * T); |
|
|
164 |
|
|
165 ## replace Ac |
|
|
166 aa(crng,crng) = matexp(crng,crng); ## discretized homegenous diff eqn |
|
|
167 |
|
|
168 ## replace Bc |
|
|
169 bb(crng,:) = matexp(crng,n+(1:m)); |
|
|
170 |
|
|
171 ## replace Acd |
|
|
172 if(nz) |
|
|
173 aa(crng,drng) = matexp(crng,n+m+(1:nz)); |
|
|
174 end |
|
|
175 |
|
|
176 stnames = sysgetsignals(sys,"st"); ## continuous states renamed below |
|
|
177 innames = sysgetsignals(sys,"in"); |
|
|
178 outnames = sysgetsignals(sys,"out"); |
|
|
179 outlist = 1:p; |
|
4771
|
180 dsys = ss(aa,bb,cc,dd,T,0,n+nz,stnames,innames, ... |
|
3431
|
181 outnames,outlist); |
|
|
182 ## rename states |
|
|
183 for ii=1:n |
|
|
184 strval = sprintf("%s_d",sysgetsignals(dsys,"st",ii,1)); |
|
|
185 dsys = syssetsignals(dsys,"st",strval,ii); |
|
|
186 endfor |
|
|
187 |
|
|
188 elseif(strcmp(opt,"bi")) |
|
|
189 if(is_digital(sys)) |
|
|
190 error("c2d: system is already digital") |
|
|
191 else |
|
|
192 ## convert with bilinear transform |
|
|
193 [a,b,c,d,tsam,n,nz,stname,inname,outname,yd] = sys2ss(sys); |
|
|
194 IT = (2/T)*eye(size(a)); |
|
|
195 A = (IT+a)/(IT-a); |
|
|
196 iab = (IT-a)\b; |
|
|
197 tk=2/sqrt(T); |
|
|
198 B = tk*iab; |
|
|
199 C = tk*(c/(IT-a)); |
|
|
200 D = d + (c*iab); |
|
|
201 stnamed = strappend(stname,"_d"); |
|
4771
|
202 dsys = ss(A,B,C,D,T,0,rows(A),stnamed,inname,outname); |
|
3775
|
203 endif |
|
3759
|
204 elseif(strcmp(opt,"matched")) |
|
|
205 if(is_digital(sys)) |
|
|
206 error("c2d: system is already digital"); |
|
|
207 elseif((length(sys.inname) != 1) || (length(sys.outname) != 1)) |
|
|
208 error("c2d: system in not single input, single output"); |
|
|
209 else |
|
|
210 sys = sysupdate(sys,"zp"); |
|
|
211 p = exp(sys.pol*T); |
|
|
212 z = exp(sys.zer*T); |
|
|
213 infinite_zeros = max(size(sys.pol))-max(size(sys.zer))-1; |
|
|
214 for i = 1:infinite_zeros |
|
|
215 z = [z ; -1]; |
|
|
216 endfor |
|
|
217 ## Should the freaquency we adjust around always be 1? |
|
|
218 [cmag,cphase,cw] = bode(sys,1); |
|
4771
|
219 [dmag,dpahse,dw] = bode(zp(z,p,1,T),1); |
|
|
220 dsys = zp(z,p,cmag/dmag,T); |
|
3431
|
221 endif |
|
|
222 else |
|
3759
|
223 error ("invalid option = %s", opt); |
|
3431
|
224 endif |
|
|
225 |
|
|
226 endfunction |
