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