Mercurial > octave-nkf

comparison scripts/control/base/__stepimp__.m @ 6448:2110cc251779

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-03-24 02:47:36 by jwe]
author jwe
date Sat, 24 Mar 2007 02:47:36 +0000
parents ddfe04062467
children 6bbf56a9718a
comparison
equal deleted inserted replaced
6447:3f79532415b5 6448:2110cc251779
35 ## Created: October 2, 1997 35 ## Created: October 2, 1997
36 ## based on lsim.m of Scottedward Hodel 36 ## based on lsim.m of Scottedward Hodel
37 37
38 function [y, t] = __stepimp__ (sitype, sys, inp, tstop, n) 38 function [y, t] = __stepimp__ (sitype, sys, inp, tstop, n)
39 39
40 if (sitype == 1) IMPULSE = 0; 40 if (sitype == 1)
41 elseif (sitype == 2) IMPULSE = 1; 41 IMPULSE = 0;
42 else error("__stepimp__: invalid sitype argument.") 42 elseif (sitype == 2)
43 endif 43 IMPULSE = 1;
44 sys = sysupdate(sys,"ss"); 44 else
45 error ("__stepimp__: invalid sitype argument");
46 endif
47 sys = sysupdate (sys, "ss");
45 48
46 USE_DEF = 0; # default tstop and n if we have to give up 49 USE_DEF = 0; # default tstop and n if we have to give up
47 N_MIN = 50; # minimum number of points 50 N_MIN = 50; # minimum number of points
48 N_MAX = 2000; # maximum number of points 51 N_MAX = 2000; # maximum number of points
49 T_DEF = 10.0; # default simulation time 52 T_DEF = 10.0; # default simulation time
50 53
51 ## collect useful information about the system 54 ## collect useful information about the system
52 [ncstates,ndstates,NIN,NOUT] = sysdimensions(sys); 55 [ncstates, ndstates, NIN, NOUT] = sysdimensions (sys);
53 TSAMPLE = sysgettsam(sys); 56 TSAMPLE = sysgettsam (sys);
54 57
55 if (nargin < 3) inp = 1; 58 if (nargin < 3)
56 elseif (inp < 1 | inp > NIN) error("Argument inp out of range") 59 inp = 1;
57 endif 60 elseif (inp < 1 || inp > NIN)
58 61 error ("__stepimp__: argument inp out of range");
59 DIGITAL = is_digital(sys); 62 endif
63
64 DIGITAL = is_digital (sys);
60 if (DIGITAL) 65 if (DIGITAL)
61 NSTATES = ndstates; 66 NSTATES = ndstates;
62 if (TSAMPLE < eps) 67 if (TSAMPLE < eps)
63 error("__stepimp__: sampling time of discrete system too small.") 68 error ("__stepimp__: sampling time of discrete system too small")
64 endif 69 endif
65 else NSTATES = ncstates; endif 70 else
71 NSTATES = ncstates;
72 endif
66 if (NSTATES < 1) 73 if (NSTATES < 1)
67 error("step: pure gain block (n_states < 1), step response is trivial"); 74 error ("__stepimp__: pure gain block (n_states < 1), step response is trivial");
68 endif 75 endif
69 if (nargin < 5) 76 if (nargin < 5)
70 ## we have to compute the time when the system reaches steady state 77 ## we have to compute the time when the system reaches steady state
71 ## and the step size 78 ## and the step size
72 ev = eig(sys2ss(sys)); 79 ev = eig (sys2ss (sys));
73 if (DIGITAL) 80 if (DIGITAL)
74 ## perform bilinear transformation on poles in z 81 ## perform bilinear transformation on poles in z
75 for i = 1:NSTATES 82 for i = 1:NSTATES
76 pole = ev(i); 83 pole = ev(i);
77 if (abs(pole + 1) < 1.0e-10) 84 if (abs(pole + 1) < 1.0e-10)
82 endfor 89 endfor
83 endif 90 endif
84 ## remove poles near zero from eigenvalue array ev 91 ## remove poles near zero from eigenvalue array ev
85 nk = NSTATES; 92 nk = NSTATES;
86 for i = 1:NSTATES 93 for i = 1:NSTATES
87 if (abs(real(ev(i))) < 1.0e-10) 94 if (abs (real (ev(i))) < 1.0e-10)
88 ev(i) = 0; 95 ev(i) = 0;
89 nk = nk - 1; 96 nk = nk - 1;
90 endif 97 endif
91 endfor 98 endfor
92 if (nk == 0) 99 if (nk == 0)
93 USE_DEF = 1; 100 USE_DEF = 1;
94 ## printf("##STEPIMP-DEBUG: using defaults.\n"); 101 ## printf("##STEPIMP-DEBUG: using defaults.\n");
95 else 102 else
96 ev = ev(find(ev)); 103 ev = ev(find (ev));
97 x = max(abs(ev)); 104 x = max (abs (ev));
98 t_step = 0.2 * pi / x; 105 t_step = 0.2 * pi / x;
99 x = min(abs(real(ev))); 106 x = min (abs (real (ev)));
100 t_sim = 5.0 / x; 107 t_sim = 5.0 / x;
101 ## round up 108 ## round up
102 yy = 10^(ceil(log10(t_sim)) - 1); 109 yy = 10^(ceil (log10 (t_sim)) - 1);
103 t_sim = yy * ceil(t_sim / yy); 110 t_sim = yy * ceil (t_sim / yy);
104 ## printf("##STEPIMP-DEBUG: nk=%d t_step=%f t_sim=%f\n", 111 ## printf("##STEPIMP-DEBUG: nk=%d t_step=%f t_sim=%f\n",
105 ## nk, t_step, t_sim); 112 ## nk, t_step, t_sim);
106 endif 113 endif
107 endif 114 endif
108 115
109 if (DIGITAL) 116 if (DIGITAL)
110 ## ---- sampled system 117 ## ---- sampled system
111 if (nargin == 5) 118 if (nargin == 5)
112 n = round(n); 119 n = round (n);
113 if (n < 2) 120 if (n < 2)
114 error("__stepimp__: n must not be less than 2.") 121 error ("__stepimp__: n must not be less than 2.")
115 endif 122 endif
116 else 123 else
117 if (nargin == 4) 124 if (nargin == 4)
118 ## n is unknown 125 ## n is unknown
119 elseif (nargin >= 1) 126 elseif (nargin >= 1)
122 tstop = (N_MIN - 1) * TSAMPLE; 129 tstop = (N_MIN - 1) * TSAMPLE;
123 else 130 else
124 tstop = t_sim; 131 tstop = t_sim;
125 endif 132 endif
126 endif 133 endif
127 n = floor(tstop / TSAMPLE) + 1; 134 n = floor (tstop / TSAMPLE) + 1;
128 if (n < 2) n = 2; endif 135 if (n < 2)
136 n = 2;
137 endif
129 if (n > N_MAX) 138 if (n > N_MAX)
130 n = N_MAX; 139 n = N_MAX;
131 printf("Hint: number of samples limited to %d by default.\n", \ 140 printf ("Hint: number of samples limited to %d by default.\n", \
132 N_MAX); 141 N_MAX);
133 printf(" ==> increase \"n\" parameter for longer simulations.\n"); 142 printf (" ==> increase \"n\" parameter for longer simulations.\n");
134 endif 143 endif
135 endif 144 endif
136 tstop = (n - 1) * TSAMPLE; 145 tstop = (n - 1) * TSAMPLE;
137 t_step = TSAMPLE; 146 t_step = TSAMPLE;
138 else 147 else
139 ## ---- continuous system 148 ## ---- continuous system
140 if (nargin == 5) 149 if (nargin == 5)
141 n = round(n); 150 n = round (n);
142 if (n < 2) 151 if (n < 2)
143 error("step: n must not be less than 2.") 152 error("step: n must not be less than 2.")
144 endif 153 endif
145 t_step = tstop / (n - 1); 154 t_step = tstop / (n - 1);
146 else 155 else
148 ## only n in unknown 157 ## only n in unknown
149 if (USE_DEF) 158 if (USE_DEF)
150 n = N_MIN; 159 n = N_MIN;
151 t_step = tstop / (n - 1); 160 t_step = tstop / (n - 1);
152 else 161 else
153 n = floor(tstop / t_step) + 1; 162 n = floor (tstop / t_step) + 1;
154 endif 163 endif
155 else 164 else
156 ## tstop and n are unknown 165 ## tstop and n are unknown
157 if (USE_DEF) 166 if (USE_DEF)
158 tstop = T_DEF; 167 tstop = T_DEF;
159 n = N_MIN; 168 n = N_MIN;
160 t_step = tstop / (n - 1); 169 t_step = tstop / (n - 1);
161 else 170 else
162 tstop = t_sim; 171 tstop = t_sim;
163 n = floor(tstop / t_step) + 1; 172 n = floor (tstop / t_step) + 1;
164 endif 173 endif
165 endif 174 endif
166 if (n < N_MIN) 175 if (n < N_MIN)
167 n = N_MIN; 176 n = N_MIN;
168 t_step = tstop / (n - 1); 177 t_step = tstop / (n - 1);
172 t_step = tstop / (N_MAX - 1); 181 t_step = tstop / (N_MAX - 1);
173 n = N_MAX; 182 n = N_MAX;
174 endif 183 endif
175 endif 184 endif
176 tstop = (n - 1) * t_step; 185 tstop = (n - 1) * t_step;
177 [jnk,B] = sys2ss(sys); 186 [jnk,B] = sys2ss (sys);
178 B = B(:,inp); 187 B = B(:,inp);
179 sys = c2d(sys, t_step); 188 sys = c2d (sys, t_step);
180 endif 189 endif
181 ## printf("##STEPIMP-DEBUG: t_step=%f n=%d tstop=%f\n", t_step, n, tstop); 190 ## printf("##STEPIMP-DEBUG: t_step=%f n=%d tstop=%f\n", t_step, n, tstop);
182 191
183 F = sys.a; 192 F = sys.a;
184 G = sys.b(:,inp); 193 G = sys.b(:,inp);
185 C = sys.c; 194 C = sys.c;
186 D = sys.d(:,inp); 195 D = sys.d(:,inp);
187 y = zeros(NOUT, n); 196 y = zeros (NOUT, n);
188 t = linspace(0, tstop, n); 197 t = linspace (0, tstop, n);
189 198
190 if (IMPULSE) 199 if (IMPULSE)
191 if (!DIGITAL && (D'*D > 0)) 200 if (! DIGITAL && D'*D > 0)
192 error("impulse: D matrix is nonzero, impulse response infinite.") 201 error ("impulse: D matrix is nonzero, impulse response infinite.")
193 endif 202 endif
194 if (DIGITAL) 203 if (DIGITAL)
195 y(:,1) = D / t_step; 204 y(:,1) = D / t_step;
196 x = G / t_step; 205 x = G / t_step;
197 else 206 else
205 endfor 214 endfor
206 if (DIGITAL) 215 if (DIGITAL)
207 y *= t_step; 216 y *= t_step;
208 endif 217 endif
209 else 218 else
210 x = zeros(NSTATES, 1); 219 x = zeros (NSTATES, 1);
211 for i = 1:n 220 for i = 1:n
212 y(:,i) = C * x + D; 221 y(:,i) = C * x + D;
213 x = F * x + G; 222 x = F * x + G;
214 endfor 223 endfor
215 endif 224 endif
216 225
217 save_automatic_replot = automatic_replot; 226 if (nargout == 0)
218 unwind_protect 227 if (IMPULSE)
219 automatic_replot(0); 228 gm = zeros (NOUT, 1);
220 if(nargout == 0) 229 tt = "impulse";
221 ## Plot the information 230 else
222 oneplot(); 231 ssys = ss (F, G, C, D, t_step);
223 __gnuplot_set__ nogrid 232 gm = dcgain (ssys);
224 __gnuplot_set__ nologscale 233 tt = "step";
225 __gnuplot_set__ autoscale 234 endif
226 __gnuplot_set__ nokey 235 ncols = floor (sqrt (NOUT));
227 if (IMPULSE) 236 nrows = ceil (NOUT / ncols);
228 gm = zeros(NOUT, 1); 237 for i = 1:NOUT
229 tt = "impulse"; 238 subplot (nrows, ncols, i);
239 if (DIGITAL)
240 [ts, ys] = stairs (t, y(i,:));
241 ts = ts(1:2*n-2)';
242 ys = ys(1:2*n-2)';
243 if (length (gm) > 0)
244 yy = [ys; gm(i)*ones(size(ts))];
245 else
246 yy = ys;
247 endif
248 plot (ts, yy);
249 grid ("on");
250 xlabel ("time [s]");
251 ylabel ("y(t)");
230 else 252 else
231 ssys = ss(F, G, C, D, t_step); 253 if (length (gm) > 0)
232 gm = dcgain(ssys); 254 yy = [y(i,:); gm(i)*ones(size(t))];
233 tt = "step";
234 endif
235 ncols = floor(sqrt(NOUT));
236 nrows = ceil(NOUT / ncols);
237 if (ncols > 1 || nrows > 1)
238 clearplot();
239 endif
240 for i = 1:NOUT
241 subplot(nrows, ncols, i);
242 title(sprintf("%s: | %s -> %s", tt,sysgetsignals(sys,"in",inp,1), ...
243 sysgetsignals(sys,"out",i,1)));
244 if (DIGITAL)
245 [ts, ys] = stairs(t, y(i,:));
246 ts = ts(1:2*n-2)'; ys = ys(1:2*n-2)';
247 if (length(gm) > 0)
248 yy = [ys; gm(i)*ones(size(ts))];
249 else
250 yy = ys;
251 endif
252 grid("on");
253 xlabel("time [s]");
254 ylabel("y(t)");
255 plot(ts, yy);
256 else 255 else
257 if (length(gm) > 0) 256 yy = y(i,:);
258 yy = [y(i,:); gm(i)*ones(size(t))];
259 else
260 yy = y(i,:);
261 endif
262 grid("on");
263 xlabel("time [s]");
264 ylabel("y(t)");
265 plot(t, yy);
266 endif 257 endif
267 endfor 258 plot (t, yy);
268 ## leave gnuplot in multiplot mode is bad style 259 grid ("on");
269 oneplot(); 260 xlabel ("time [s]");
270 y=[]; 261 ylabel ("y(t)");
271 t=[]; 262 endif
272 endif 263 title (sprintf ("%s: | %s -> %s", tt,
273 ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n"); 264 sysgetsignals (sys, "in", inp, 1),
274 unwind_protect_cleanup 265 sysgetsignals (sys, "out", i, 1)));
275 automatic_replot(save_automatic_replot); 266 endfor
276 end_unwind_protect 267 y = [];
268 t = [];
269 endif
270 ## printf("##STEPIMP-DEBUG: gratulations, successfull completion.\n");
277 endfunction 271 endfunction