octave_packages/control-2.3.52/__time_response__.m

   1 ## Copyright (C) 2009, 2010   Lukas F. Reichlin
   2 ##
   3 ## This file is part of LTI Syncope.
   4 ##
   5 ## LTI Syncope is free software: you can redistribute it and/or modify
   6 ## it under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation, either version 3 of the License, or
   8 ## (at your option) any later version.
   9 ##
  10 ## LTI Syncope is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
  17
  18 ## -*- texinfo -*-
  19 ## Common code for the time response functions step, impulse and initial.
  20
  21 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  22 ## Created: October 2009
  23 ## Version: 0.2
  24
  25 function [y, t, x_arr] = __time_response__ (sys, resptype, plotflag, tfinal, dt, x0, sysname)
  26
  27   if (! isa (sys, "ss"))
  28     sys = ss (sys);                                    # sys must be proper
  29   endif
  30
  31   if (is_real_vector (tfinal) && length (tfinal) > 1)  # time vector t passed
  32     dt = tfinal(2) - tfinal(1);                        # assume that t is regularly spaced
  33     tfinal = tfinal(end);
  34   endif
  35
  36   [A, B, C, D, tsam] = ssdata (sys);
  37
  38   discrete = ! isct (sys);                             # static gains are treated as analog systems
  39   tsam = abs (tsam);                                   # use 1 second if tsam is unspecified (-1)
  40
  41   if (discrete)
  42     if (! isempty (dt))
  43       warning ("time_response: argument dt has no effect on sampling time of discrete system");
  44     endif
  45
  46     dt = tsam;
  47   endif
  48
  49   [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, dt);
  50
  51   if (! discrete)
  52     sys = c2d (sys, dt, "zoh");
  53   endif
  54
  55   [F, G] = ssdata (sys);                               # matrices C and D don't change
  56
  57   n = rows (F);                                        # number of states
  58   m = columns (G);                                     # number of inputs
  59   p = rows (C);                                        # number of outputs
  60
  61   ## time vector
  62   t = reshape (0 : dt : tfinal, [], 1);
  63   l_t = length (t);
  64
  65   switch (resptype)
  66     case "initial"
  67       str = ["Response of ", sysname, " to Initial Conditions"];
  68       yfinal = zeros (p, 1);
  69
  70       ## preallocate memory
  71       y = zeros (l_t, p);
  72       x_arr = zeros (l_t, n);
  73
  74       ## initial conditions
  75       x = reshape (x0, [], 1);                         # make sure that x is a column vector
  76
  77       if (n != length (x0) || ! is_real_vector (x0))
  78         error ("initial: x0 must be a real vector with %d elements", n);
  79       endif
  80
  81       ## simulation
  82       for k = 1 : l_t
  83         y(k, :) = C * x;
  84         x_arr(k, :) = x;
  85         x = F * x;
  86       endfor
  87
  88     case "step"
  89       str = ["Step Response of ", sysname];
  90       yfinal = dcgain (sys);
  91
  92       ## preallocate memory
  93       y = zeros (l_t, p, m);
  94       x_arr = zeros (l_t, n, m);
  95
  96       for j = 1 : m                                    # for every input channel
  97         ## initial conditions
  98         x = zeros (n, 1);
  99         u = zeros (m, 1);
 100         u(j) = 1;
 101
 102         ## simulation
 103         for k = 1 : l_t
 104           y(k, :, j) = C * x + D * u;
 105           x_arr(k, :, j) = x;
 106           x = F * x + G * u;
 107         endfor
 108       endfor
 109
 110     case "impulse"
 111       str = ["Impulse Response of ", sysname];
 112       yfinal = zeros (p, m);
 113
 114       ## preallocate memory
 115       y = zeros (l_t, p, m);
 116       x_arr = zeros (l_t, n, m);
 117
 118       for j = 1 : m                                    # for every input channel
 119         ## initial conditions
 120         u = zeros (m, 1);
 121         u(j) = 1;
 122
 123         if (discrete)
 124           x = zeros (n, 1);                            # zero by definition
 125           y(1, :, j) = D * u / dt;
 126           x_arr(1, :, j) = x;
 127           x = G * u / dt;
 128         else
 129           x = B * u;                                   # B, not G!
 130           y(1, :, j) = C * x;
 131           x_arr(1, :, j) = x;
 132           x = F * x;
 133         endif
 134
 135         ## simulation
 136         for k = 2 : l_t
 137           y (k, :, j) = C * x;
 138           x_arr(k, :, j) = x;
 139           x = F * x;
 140         endfor
 141       endfor
 142
 143       if (discrete)
 144         y *= dt;
 145         x_arr *= dt;
 146       endif
 147
 148     otherwise
 149       error ("time_response: invalid response type");
 150
 151   endswitch
 152
 153
 154   if (plotflag)                                        # display plot
 155
 156     ## TODO: Set correct titles, especially for multi-input systems
 157
 158     stable = isstable (sys);
 159     outname = get (sys, "outname");
 160     outname = __labels__ (outname, "y_");
 161
 162     if (strcmp (resptype, "initial"))
 163       cols = 1;
 164     else
 165       cols = m;
 166     endif
 167
 168     if (discrete)                                      # discrete system
 169       for k = 1 : p
 170         for j = 1 : cols
 171
 172           subplot (p, cols, (k-1)*cols+j);
 173
 174           if (stable)
 175             stairs (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
 176           else
 177             stairs (t, y(:, k, j));
 178           endif
 179
 180           grid ("on");
 181
 182           if (k == 1 && j == 1)
 183             title (str);
 184           endif
 185
 186           if (j == 1)
 187             ylabel (sprintf ("Amplitude %s", outname{k}));
 188           endif
 189
 190         endfor
 191       endfor
 192
 193       xlabel ("Time [s]");
 194
 195     else                                               # continuous system
 196       for k = 1 : p
 197         for j = 1 : cols
 198
 199           subplot (p, cols, (k-1)*cols+j);
 200
 201           if (stable)
 202             plot (t, [y(:, k, j), yfinal(k, j) * ones(l_t, 1)]);
 203           else
 204             plot (t, y(:, k, j));
 205           endif
 206
 207           grid ("on");
 208
 209           if (k == 1 && j == 1)
 210             title (str);
 211           endif
 212
 213           if (j == 1)
 214             ylabel (sprintf ("Amplitude %s", outname{k}));
 215           endif
 216
 217         endfor
 218       endfor
 219
 220       xlabel ("Time [s]");
 221
 222     endif
 223   endif
 224
 225 endfunction
 226
 227
 228 function [tfinal, dt] = __sim_horizon__ (A, discrete, tfinal, Ts)
 229
 230   ## code based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel
 231
 232   TOL = 1.0e-10;                                       # values below TOL are assumed to be zero
 233   N_MIN = 50;                                          # min number of points
 234   N_MAX = 2000;                                        # max number of points
 235   N_DEF = 1000;                                        # default number of points
 236   T_DEF = 10;                                          # default simulation time
 237
 238   n = rows (A);
 239   eigw = eig (A);
 240
 241   if (discrete)
 242     ## perform bilinear transformation on poles in z
 243     for k = 1 : n
 244       pol = eigw(k);
 245       if (abs (pol + 1) < TOL)
 246         eigw(k) = 0;
 247       else
 248         eigw(k) = 2 / Ts * (pol - 1) / (pol + 1);
 249       endif
 250     endfor
 251   endif
 252
 253   ## remove poles near zero from eigenvalue array eigw
 254   nk = n;
 255   for k = 1 : n
 256     if (abs (real (eigw(k))) < TOL)
 257       eigw(k) = 0;
 258       nk -= 1;
 259     endif
 260   endfor
 261
 262   if (nk == 0)
 263     if (isempty (tfinal))
 264       tfinal = T_DEF;
 265     endif
 266
 267     if (! discrete)
 268       dt = tfinal / N_DEF;
 269     endif
 270   else
 271     eigw = eigw(find (eigw));
 272     eigw_max = max (abs (eigw));
 273
 274     if (! discrete)
 275       dt = 0.2 * pi / eigw_max;
 276     endif
 277
 278     if (isempty (tfinal))
 279       eigw_min = min (abs (real (eigw)));
 280       tfinal = 5.0 / eigw_min;
 281
 282       ## round up
 283       yy = 10^(ceil (log10 (tfinal)) - 1);
 284       tfinal = yy * ceil (tfinal / yy);
 285     endif
 286
 287     if (! discrete)
 288       N = tfinal / dt;
 289
 290       if (N < N_MIN)
 291         dt = tfinal / N_MIN;
 292       endif
 293
 294       if (N > N_MAX)
 295         dt = tfinal / N_MAX;
 296       endif
 297     endif
 298   endif
 299
 300   if (! isempty (Ts))                                  # catch case cont. system with dt specified
 301     dt = Ts;
 302   endif
 303
 304 endfunction