1 ## Copyright (C) 2009, 2010, 2011 Lukas F. Reichlin
3 ## This file is part of LTI Syncope.
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.
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.
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/>.
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.
30 ## LTI model. System must be proper, i.e. it must not have more zeros than poles.
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.
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)}.
41 ## Vector of initial conditions for each state. If not specified, a zero vector is assumed.
43 ## Discretization method for continuous-time models. Default value is zoh
44 ## (zero-order hold). All methods from @code{c2d} are supported.
50 ## Output response array. Has as many rows as time samples (length of t)
51 ## and as many columns as outputs.
53 ## Time row vector. It is always evenly spaced.
55 ## State trajectories array. Has @code{length (t)} rows and as many columns as states.
58 ## @seealso{impulse, initial, step}
61 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
62 ## Created: October 2009
65 function [y_r, t_r, x_r] = lsim (sys, u, t = [], x0 = [], method = "zoh")
67 ## TODO: multiplot feature: lsim (sys1, "b", sys2, "r", ..., u, t)
69 if (nargin < 2 || nargin > 5)
73 if (! isa (sys, "ss"))
74 sys = ss (sys); # sys must be proper
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");
83 if (! is_real_vector (t) && ! isempty (t))
84 error ("lsim: time vector t must be real or empty");
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)
92 if (discrete) # discrete system
93 if (isempty (t)) # lsim (sys, u)
96 tfinal = tsam * (urows - 1);
97 elseif (length (t) == 1) # lsim (sys, u, tfinal)
101 else # lsim (sys, u, t, ...)
102 warning ("lsim: spacing of time vector has no effect on sampling time of discrete system");
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);
114 else # lsim (sys, u, t, ...)
115 dt = t(2) - t(1); # assume that t is regularly spaced
119 sys = c2d (sys, dt, method); # convert to discrete-time model
122 [A, B, C, D] = ssdata (sys);
124 n = rows (A); # number of states
125 m = columns (B); # number of inputs
126 p = rows (C); # number of outputs
128 t = reshape (tinitial : dt : tfinal, [], 1); # time vector
132 error ("lsim: input vector u must have %d rows", trows);
136 error ("lsim: input vector u must have %d columns", m);
139 ## preallocate memory
140 y = zeros (trows, p);
141 x_arr = zeros (trows, n);
143 ## initial conditions
146 elseif (n != length (x0) || ! is_real_vector (x0))
147 error ("initial: x0 must be a vector with %d elements", n);
150 x = reshape (x0, [], 1); # make sure that x is a column vector
154 y(k, :) = C * x + D * u(k, :).';
156 x = A * x + B * u(k, :).';
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
171 ylabel (sprintf ("Amplitude %s", outname{k}));
174 else # continuous system
182 ylabel (sprintf ("Amplitude %s", outname{k}));
195 ## TODO: add test cases