octave_packages/control-2.3.52/lsim.m

   1 ## Copyright (C) 2009, 2010, 2011   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 ## @deftypefn{Function File} {[@var{y}, @var{t}, @var{x}] =} lsim (@var{sys}, @var{u})
  20 ## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} lsim (@var{sys}, @var{u}, @var{t})
  21 ## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} lsim (@var{sys}, @var{u}, @var{t}, @var{x0})
  22 ## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} lsim (@var{sys}, @var{u}, @var{t}, @var{[]}, @var{method})
  23 ## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} lsim (@var{sys}, @var{u}, @var{t}, @var{x0}, @var{method})
  24 ## Simulate LTI model response to arbitrary inputs.  If no output arguments are given,
  25 ## the system response is plotted on the screen.
  26 ##
  27 ## @strong{Inputs}
  28 ## @table @var
  29 ## @item sys
  30 ## LTI model.  System must be proper, i.e. it must not have more zeros than poles.
  31 ## @item u
  32 ## Vector or array of input signal.  Needs @code{length(t)} rows and as many columns
  33 ## as there are inputs.  If @var{sys} is a single-input system, row vectors @var{u}
  34 ## of length @code{length(t)} are accepted as well.
  35 ## @item t
  36 ## Time vector.  Should be evenly spaced.  If @var{sys} is a continuous-time system
  37 ## and @var{t} is a real scalar, @var{sys} is discretized with sampling time
  38 ## @code{tsam = t/(rows(u)-1)}.  If @var{sys} is a discrete-time system and @var{t}
  39 ## is not specified, vector @var{t} is assumed to be @code{0 : tsam : tsam*(rows(u)-1)}.
  40 ## @item x0
  41 ## Vector of initial conditions for each state.  If not specified, a zero vector is assumed.
  42 ## @item method
  43 ## Discretization method for continuous-time models.  Default value is zoh
  44 ## (zero-order hold).  All methods from @code{c2d} are supported.
  45 ## @end table
  46 ##
  47 ## @strong{Outputs}
  48 ## @table @var
  49 ## @item y
  50 ## Output response array.  Has as many rows as time samples (length of t)
  51 ## and as many columns as outputs.
  52 ## @item t
  53 ## Time row vector.  It is always evenly spaced.
  54 ## @item x
  55 ## State trajectories array.  Has @code{length (t)} rows and as many columns as states.
  56 ## @end table
  57 ##
  58 ## @seealso{impulse, initial, step}
  59 ## @end deftypefn
  60
  61 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  62 ## Created: October 2009
  63 ## Version: 0.3
  64
  65 function [y_r, t_r, x_r] = lsim (sys, u, t = [], x0 = [], method = "zoh")
  66
  67   ## TODO: multiplot feature:   lsim (sys1, "b", sys2, "r", ..., u, t)
  68
  69   if (nargin < 2 || nargin > 5)
  70     print_usage ();
  71   endif
  72
  73   if (! isa (sys, "ss"))
  74     sys = ss (sys);                             # sys must be proper
  75   endif
  76
  77   if (is_real_vector (u))
  78     u = reshape (u, [], 1);                     # allow row vectors for single-input systems
  79   elseif (! is_real_matrix (u))
  80     error ("lsim: input signal u must be an array of real numbers");
  81   endif
  82
  83   if (! is_real_vector (t) && ! isempty (t))
  84     error ("lsim: time vector t must be real or empty");
  85   endif
  86
  87   discrete = ! isct (sys);                      # static gains are treated as continuous-time systems
  88   tsam = abs (get (sys, "tsam"));               # use 1 second as default if tsam is unspecified (-1)
  89   urows = rows (u);
  90   ucols = columns (u);
  91
  92   if (discrete)                                 # discrete system
  93     if (isempty (t))                            # lsim (sys, u)
  94       dt = tsam;
  95       tinitial = 0;
  96       tfinal = tsam * (urows - 1);
  97     elseif (length (t) == 1)                    # lsim (sys, u, tfinal)
  98       dt = tsam;
  99       tinitial = 0;
 100       tfinal = t;
 101     else                                        # lsim (sys, u, t, ...)
 102       warning ("lsim: spacing of time vector has no effect on sampling time of discrete system");
 103       dt = tsam;
 104       tinitial = t(1);
 105       tfinal = t(end);
 106     endif
 107   else                                          # continuous system
 108     if (isempty (t))                            # lsim (sys, u, [], ...)
 109       error ("lsim: invalid time vector");
 110     elseif (length (t) == 1)                    # lsim (sys, u, tfinal, ...)
 111       dt = t / (urows - 1);
 112       tinitial = 0;
 113       tfinal = t;
 114     else                                        # lsim (sys, u, t, ...)
 115       dt = t(2) - t(1);                         # assume that t is regularly spaced
 116       tinitial = t(1);
 117       tfinal = t(end);
 118     endif
 119     sys = c2d (sys, dt, method);                # convert to discrete-time model
 120   endif
 121
 122   [A, B, C, D] = ssdata (sys);
 123
 124   n = rows (A);                                 # number of states
 125   m = columns (B);                              # number of inputs
 126   p = rows (C);                                 # number of outputs
 127
 128   t = reshape (tinitial : dt : tfinal, [], 1);  # time vector
 129   trows = length (t);
 130
 131   if (urows != trows)
 132     error ("lsim: input vector u must have %d rows", trows);
 133   endif
 134
 135   if (ucols != m)
 136     error ("lsim: input vector u must have %d columns", m);
 137   endif
 138
 139   ## preallocate memory
 140   y = zeros (trows, p);
 141   x_arr = zeros (trows, n);
 142
 143   ## initial conditions
 144   if (isempty (x0))
 145     x0 = zeros (n, 1);
 146   elseif (n != length (x0) || ! is_real_vector (x0))
 147     error ("initial: x0 must be a vector with %d elements", n);
 148   endif
 149
 150   x = reshape (x0, [], 1);                      # make sure that x is a column vector
 151
 152   ## simulation
 153   for k = 1 : trows
 154     y(k, :) = C * x  +  D * u(k, :).';
 155     x_arr(k, :) = x;
 156     x = A * x  +  B * u(k, :).';
 157   endfor
 158
 159   if (nargout == 0)                             # plot information
 160     str = ["Linear Simulation Results of ", inputname(1)];
 161     outname = get (sys, "outname");
 162     outname = __labels__ (outname, "y_");
 163     if (discrete)                               # discrete system
 164       for k = 1 : p
 165         subplot (p, 1, k);
 166         stairs (t, y(:, k));
 167         grid ("on");
 168         if (k == 1)
 169           title (str);
 170         endif
 171         ylabel (sprintf ("Amplitude %s", outname{k}));
 172       endfor
 173       xlabel ("Time [s]");
 174     else                                        # continuous system
 175       for k = 1 : p
 176         subplot (p, 1, k);
 177         plot (t, y(:, k));
 178         grid ("on");
 179         if (k == 1)
 180           title (str);
 181         endif
 182         ylabel (sprintf ("Amplitude %s", outname{k}));
 183       endfor
 184       xlabel ("Time [s]");
 185     endif
 186   else                                          # return values
 187     y_r = y;
 188     t_r = t;
 189     x_r = x_arr;
 190   endif
 191
 192 endfunction
 193
 194
 195 ## TODO: add test cases