octave_packages/control-2.3.52/covar.m

   1 ## Copyright (C) 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 ## @deftypefn{Function File} {[@var{p}, @var{q}] =} covar (@var{sys}, @var{w})
  20 ## Return the steady-state covariance.
  21 ##
  22 ## @strong{Inputs}
  23 ## @table @var
  24 ## @item sys
  25 ## LTI model.
  26 ## @item w
  27 ## Intensity of Gaussian white noise inputs which drive @var{sys}.
  28 ## @end table
  29 ##
  30 ## @strong{Outputs}
  31 ## @table @var
  32 ## @item p
  33 ## Output covariance.
  34 ## @item q
  35 ## State covariance.
  36 ## @end table
  37 ##
  38 ## @seealso{lyap, dlyap}
  39 ## @end deftypefn
  40
  41 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  42 ## Created: January 2010
  43 ## Version: 0.1
  44
  45 function [p, q] = covar (sys, w)
  46
  47   if (nargin != 2)
  48     print_usage ();
  49   endif
  50
  51   if (! isa (sys, "lti"))
  52     error ("covar: first argument must be an LTI model");
  53   endif
  54
  55   if (! isstable (sys))
  56     error ("covar: system must be stable");
  57   endif
  58
  59   [a, b, c, d] = ssdata (sys);
  60
  61   if (isct (sys))
  62     if (any (d(:)))
  63       error ("covar: system is not strictly proper");
  64     endif
  65
  66     q = lyap (a, b*w*b.');
  67     p = c*q*c.';
  68   else
  69     q = dlyap (a, b*w*b.');
  70     p = c*q*c.' + d*w*d.';
  71   endif
  72
  73 endfunction
  74
  75 ## continuous-time
  76 %!shared p, q, p_exp, q_exp
  77 %! sys = ss (-1, 1, 1, 0);
  78 %! [p, q] = covar (sys, 5);
  79 %! p_exp = 2.5000;
  80 %! q_exp = 2.5000;
  81 %!assert (p, p_exp, 1e-4);
  82 %!assert (q, q_exp, 1e-4);
  83
  84 ## discrete-time
  85 %!shared p, q, p_exp, q_exp
  86 %! sys = ss ([-0.2, -0.5; 1, 0], [2; 0], [1, 0.5], [0], 0.1);
  87 %! [p, q] = covar (sys, 5);
  88 %! p_exp = 30.3167;
  89 %! q_exp = [27.1493, -3.6199; -3.6199, 27.1493];
  90 %!assert (p, p_exp, 1e-4);
  91 %!assert (q, q_exp, 1e-4);