octave_packages/control-2.3.52/@lti/norm.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{gain} =} norm (@var{sys}, @var{2})
  20 ## @deftypefnx {Function File} {[@var{gain}, @var{wpeak}] =} norm (@var{sys}, @var{inf})
  21 ## @deftypefnx {Function File} {[@var{gain}, @var{wpeak}] =} norm (@var{sys}, @var{inf}, @var{tol})
  22 ## Return H-2 or L-inf norm of LTI model.
  23 ##
  24 ## @strong{Algorithm}@*
  25 ## Uses SLICOT AB13BD and AB13DD by courtesy of
  26 ## @uref{http://www.slicot.org, NICONET e.V.}
  27 ## @end deftypefn
  28
  29 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  30 ## Created: November 2009
  31 ## Version: 0.5
  32
  33 function [gain, varargout] = norm (sys, ntype = "2", tol = 0.01)
  34
  35   if (nargin > 3)                     # norm () is caught by built-in function
  36     print_usage ();
  37   endif
  38
  39   if (is_real_scalar (ntype))
  40     if (ntype == 2)
  41       ntype = "2";
  42     elseif (isinf (ntype))
  43       ntype = "inf";
  44     else
  45       error ("lti: norm: invalid norm type");
  46     endif
  47   elseif (ischar (ntype))
  48     ntype = lower (ntype);
  49   else
  50     error ("lti: norm: invalid norm type");
  51   endif
  52
  53   switch (ntype)
  54     case "2"
  55       gain = h2norm (sys);
  56
  57     case "inf"
  58       [gain, varargout{1}] = linfnorm (sys, tol);
  59
  60     otherwise
  61       error ("lti: norm: invalid norm type");
  62   endswitch
  63
  64 endfunction
  65
  66
  67 function gain = h2norm (sys)
  68
  69   if (isstable (sys))
  70     [a, b, c, d] = ssdata (sys);
  71     discrete = ! isct (sys);
  72     if (! discrete && any (d(:)))     # continuous and non-zero feedthrough
  73       gain = inf;
  74     else
  75       gain = slab13bd (a, b, c, d, discrete);
  76     endif
  77   else
  78     gain = inf;
  79   endif
  80
  81 endfunction
  82
  83
  84 function [gain, wpeak] = linfnorm (sys, tol = 0.01)
  85
  86   [a, b, c, d, e, tsam, scaled] = dssdata (sys, []);
  87   discrete = ! isct (sys);
  88   tol = max (tol, 100*eps);
  89
  90   if (isempty (e))
  91     [fpeak, gpeak] = slab13dd (a, a, b, c, d, discrete, false, tol, scaled);  # TODO: avoid dummy argument
  92   else
  93     if (rcond (e) < eps)
  94       gain = inf;
  95       wpeak = inf;
  96       return;
  97     else
  98       [fpeak, gpeak] = slab13dd (a, e, b, c, d, discrete, true, tol, scaled);
  99     endif
 100   endif
 101
 102   if (fpeak(2) > 0)
 103     if (discrete)
 104       wpeak = fpeak(1) / abs (tsam);  # tsam could be -1
 105     else
 106       wpeak = fpeak(1);
 107     endif
 108   else
 109     wpeak = inf;
 110   endif
 111
 112   if (gpeak(2) > 0)
 113     gain = gpeak(1);
 114   else
 115     gain = inf;
 116   endif
 117
 118 endfunction
 119
 120
 121 ## norm ct
 122 %!shared H2, Hinf
 123 %! sys = ss (-1, 1, 1, 0);
 124 %! H2 = norm (sys, 2);
 125 %! Hinf = norm (sys, inf);
 126 %!assert (H2, 0.7071, 1.5e-5);
 127 %!assert (Hinf, 1, 5e-4);
 128
 129
 130 ## norm dt
 131 %!shared H2, Hinf
 132 %! a = [ 2.417   -1.002    0.5488
 133 %!           2        0         0
 134 %!           0      0.5         0 ];
 135 %! b = [     1
 136 %!           0
 137 %!           0 ];
 138 %! c = [-0.424    0.436   -0.4552 ];
 139 %! d = [     1 ];
 140 %! sys = ss (a, b, c, d, 0.1);
 141 %! H2 = norm (sys, 2);
 142 %! Hinf = norm (sys, inf);
 143 %!assert (H2, 1.2527, 1.5e-5);
 144 %!assert (Hinf, 2.7, 0.1);