]> Creatis software - CreaPhase.git/blobdiff - octave_packages/control-2.3.52/initial.m
Add a useful package (from Source forge) for octave
[CreaPhase.git] / octave_packages / control-2.3.52 / initial.m
diff --git a/octave_packages/control-2.3.52/initial.m b/octave_packages/control-2.3.52/initial.m
new file mode 100644 (file)
index 0000000..b744813
--- /dev/null
@@ -0,0 +1,151 @@
+## Copyright (C) 2009   Lukas F. Reichlin
+##
+## This file is part of LTI Syncope.
+##
+## LTI Syncope is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## LTI Syncope is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with LTI Syncope.  If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0})
+## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}, @var{t})
+## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}, @var{tfinal})
+## @deftypefnx{Function File} {[@var{y}, @var{t}, @var{x}] =} initial (@var{sys}, @var{x0}, @var{tfinal}, @var{dt})
+## Initial condition response of state-space model.
+## If no output arguments are given, the response is printed on the screen.
+##
+## @strong{Inputs}
+## @table @var
+## @item sys
+## State-space model.
+## @item x0
+## Vector of initial conditions for each state.
+## @item t
+## Optional time vector.  Should be evenly spaced.  If not specified, it is calculated
+## by the poles of the system to reflect adequately the response transients.
+## @item tfinal
+## Optional simulation horizon.  If not specified, it is calculated by
+## the poles of the system to reflect adequately the response transients.
+## @item dt
+## Optional sampling time.  Be sure to choose it small enough to capture transient
+## phenomena.  If not specified, it is calculated by the poles of the system.
+## @end table
+##
+## @strong{Outputs}
+## @table @var
+## @item y
+## Output response array.  Has as many rows as time samples (length of t)
+## and as many columns as outputs.
+## @item t
+## Time row vector.
+## @item x
+## State trajectories array.  Has @code{length (t)} rows and as many columns as states.
+## @end table
+##
+## @strong{Example}
+## @example
+## @group
+##                    .
+## Continuous Time:   x = A x ,   y = C x ,   x(0) = x0
+##
+## Discrete Time:   x[k+1] = A x[k] ,   y[k] = C x[k] ,   x[0] = x0
+## @end group
+## @end example
+##
+## @seealso{impulse, lsim, step}
+## @end deftypefn
+
+## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
+## Created: October 2009
+## Version: 0.1
+
+function [y_r, t_r, x_r] = initial (sys, x0, tfinal = [], dt = [])
+
+  ## TODO: multiplot feature:   initial (sys1, "b", sys2, "r", ..., x0, ...)
+
+  if (nargin < 2 || nargin > 4)
+    print_usage ();
+  endif
+
+  [y, t, x] = __time_response__ (sys, "initial", ! nargout, tfinal, dt, x0, inputname (1));
+
+  if (nargout)
+    y_r = y;
+    t_r = t;
+    x_r = x;
+  endif
+
+endfunction
+
+
+%!shared initial_c, initial_c_exp, initial_d, initial_d_exp
+%!
+%! A = [ -2.8    2.0   -1.8
+%!       -2.4   -2.0    0.8
+%!        1.1    1.7   -1.0 ];
+%!
+%! B = [ -0.8    0.5    0
+%!        0      0.7    2.3
+%!       -0.3   -0.1    0.5 ];
+%!
+%! C = [ -0.1    0     -0.3
+%!        0.9    0.5    1.2
+%!        0.1   -0.1    1.9 ];
+%!
+%! D = [ -0.5    0      0
+%!        0.1    0      0.3
+%!       -0.8    0      0   ];
+%!
+%! x_0 = [1, 2, 3];
+%!
+%! sysc = ss (A, B, C, D);
+%!
+%! [yc, tc, xc] = initial (sysc, x_0, 0.2, 0.1);
+%! initial_c = [yc, tc, xc];
+%!
+%! sysd = c2d (sysc, 2);
+%!
+%! [yd, td, xd] = initial (sysd, x_0, 4);
+%! initial_d = [yd, td, xd];
+%!
+%! ## expected values computed by the "dark side"
+%!
+%! yc_exp = [ -1.0000    5.5000    5.6000
+%!            -0.9872    5.0898    5.7671
+%!            -0.9536    4.6931    5.7598 ];
+%!
+%! tc_exp = [  0.0000
+%!             0.1000
+%!             0.2000 ];
+%!
+%! xc_exp = [  1.0000    2.0000    3.0000
+%!             0.5937    1.6879    3.0929
+%!             0.2390    1.5187    3.0988 ];
+%!
+%! initial_c_exp = [yc_exp, tc_exp, xc_exp];
+%!
+%! yd_exp = [ -1.0000    5.5000    5.6000
+%!            -0.6550    3.1673    4.2228
+%!            -0.5421    2.6186    3.4968 ];
+%!
+%! td_exp = [  0
+%!             2
+%!             4 ];
+%!
+%! xd_exp = [  1.0000    2.0000    3.0000
+%!            -0.4247    1.5194    2.3249
+%!            -0.3538    1.2540    1.9250 ];
+%!
+%! initial_d_exp = [yd_exp, td_exp, xd_exp];
+%!
+%!assert (initial_c, initial_c_exp, 1e-4)
+%!assert (initial_d, initial_d_exp, 1e-4)