octave_packages/control-2.3.52/lyap.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{x} =} lyap (@var{a}, @var{b})
  20 ## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{c})
  21 ## @deftypefnx{Function File} {@var{x} =} lyap (@var{a}, @var{b}, @var{[]}, @var{e})
  22 ## Solve continuous-time Lyapunov or Sylvester equations.
  23 ##
  24 ## @strong{Equations}
  25 ## @example
  26 ## @group
  27 ## AX + XA' + B = 0      (Lyapunov Equation)
  28 ##
  29 ## AX + XB + C = 0       (Sylvester Equation)
  30 ##
  31 ## AXE' + EXA' + B = 0   (Generalized Lyapunov Equation)
  32 ## @end group
  33 ## @end example
  34 ##
  35 ## @strong{Algorithm}@*
  36 ## Uses SLICOT SB03MD, SB04MD and SG03AD by courtesy of
  37 ## @uref{http://www.slicot.org, NICONET e.V.}
  38 ##
  39 ## @seealso{lyapchol, dlyap, dlyapchol}
  40 ## @end deftypefn
  41
  42 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  43 ## Created: January 2010
  44 ## Version: 0.2.1
  45
  46 function [x, scale] = lyap (a, b, c, e)
  47
  48   scale = 1;
  49
  50   switch (nargin)
  51     case 2                                      # Lyapunov equation
  52
  53       if (! is_real_square_matrix (a, b))
  54         ## error ("lyap: a, b must be real and square");
  55         error ("lyap: %s, %s must be real and square", \
  56                 inputname (1), inputname (2));
  57       endif
  58
  59       if (rows (a) != rows (b))
  60         ## error ("lyap: a, b must have the same number of rows");
  61         error ("lyap: %s, %s must have the same number of rows", \
  62                 inputname (1), inputname (2));
  63
  64       endif
  65
  66       [x, scale] = slsb03md (a, -b, false);     # AX + XA' = -B
  67
  68       ## x /= scale;                            # 0 < scale <= 1
  69
  70     case 3                                      # Sylvester equation
  71
  72       if (! is_real_square_matrix (a, b))
  73         ## error ("lyap: a, b must be real and square");
  74         error ("lyap: %s, %s must be real and square", \
  75                 inputname (1), inputname (2));
  76       endif
  77
  78       if (! is_real_matrix (c) || rows (c) != rows (a) || columns (c) != columns (b))
  79         ## error ("lyap: c must be a real (%dx%d) matrix", rows (a), columns (b));
  80         error ("lyap: %s must be a real (%dx%d) matrix", \
  81                 rows (a), columns (b), inputname (3));
  82       endif
  83
  84       x = slsb04md (a, b, -c);  # AX + XB = -C
  85
  86     case 4                                      # generalized Lyapunov equation
  87
  88       if (! isempty (c))
  89         print_usage ();
  90       endif
  91
  92       if (! is_real_square_matrix (a, b, e))
  93         ## error ("lyap: a, b, e must be real and square");
  94         error ("lyap: %s, %s, %s must be real and square", \
  95                 inputname (1), inputname (2), inputname (4));
  96       endif
  97
  98       if (rows (b) != rows (a) || rows (e) != rows (a))
  99         ## error ("lyap: a, b, e must have the same number of rows");
 100         error ("lyap: %s, %s, %s must have the same number of rows", \
 101                 inputname (1), inputname (2), inputname (4));
 102       endif
 103
 104       if (! issymmetric (b))
 105         ## error ("lyap: b must be symmetric");
 106         error ("lyap: %s must be symmetric", \
 107                 inputname (2));
 108       endif
 109
 110       [x, scale] = slsg03ad (a, e, -b, false);  # AXE' + EXA' = -B
 111
 112       ## x /= scale;                            # 0 < scale <= 1
 113
 114     otherwise
 115       print_usage ();
 116
 117   endswitch
 118
 119   if (scale < 1)
 120     warning ("lyap: solution scaled by %g to prevent overflow", scale);
 121   endif
 122
 123 endfunction
 124
 125
 126 ## Lyapunov
 127 %!shared X, X_exp
 128 %! A = [1, 2; -3, -4];
 129 %! Q = [3, 1; 1, 1];
 130 %! X = lyap (A, Q);
 131 %! X_exp = [ 6.1667, -3.8333;
 132 %!          -3.8333,  3.0000];
 133 %!assert (X, X_exp, 1e-4);
 134
 135 ## Sylvester
 136 %!shared X, X_exp
 137 %! A = [2.0   1.0   3.0
 138 %!      0.0   2.0   1.0
 139 %!      6.0   1.0   2.0];
 140 %!
 141 %! B = [2.0   1.0
 142 %!      1.0   6.0];
 143 %!
 144 %! C = [2.0   1.0
 145 %!      1.0   4.0
 146 %!      0.0   5.0];
 147 %!
 148 %! X = lyap (A, B, -C);
 149 %!
 150 %! X_exp = [-2.7685   0.5498
 151 %!          -1.0531   0.6865
 152 %!           4.5257  -0.4389];
 153 %!
 154 %!assert (X, X_exp, 1e-4);
 155
 156 ## Generalized Lyapunov
 157 %!shared X, X_exp
 158 %! A = [  3.0     1.0     1.0
 159 %!        1.0     3.0     0.0
 160 %!        1.0     0.0     2.0];
 161 %!
 162 %! E = [  1.0     3.0     0.0
 163 %!        3.0     2.0     1.0
 164 %!        1.0     0.0     1.0];
 165 %!
 166 %! B = [-64.0   -73.0   -28.0
 167 %!      -73.0   -70.0   -25.0
 168 %!      -28.0   -25.0   -18.0];
 169 %!
 170 %! X = lyap (A.', -B, [], E.');
 171 %!
 172 %! X_exp = [-2.0000  -1.0000   0.0000
 173 %!          -1.0000  -3.0000  -1.0000
 174 %!           0.0000  -1.0000  -3.0000];
 175 %!
 176 %!assert (X, X_exp, 1e-4);
 177