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