octave_packages/control-2.3.52/dlyap.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} =} dlyap (@var{a}, @var{b})
  20 ## @deftypefnx{Function File} {@var{x} =} dlyap (@var{a}, @var{b}, @var{c})
  21 ## @deftypefnx{Function File} {@var{x} =} dlyap (@var{a}, @var{b}, @var{[]}, @var{e})
  22 ## Solve discrete-time Lyapunov or Sylvester equations.
  23 ##
  24 ## @strong{Equations}
  25 ## @example
  26 ## @group
  27 ## AXA' - X + B = 0      (Lyapunov Equation)
  28 ##
  29 ## AXB' - X + C = 0      (Sylvester Equation)
  30 ##
  31 ## AXA' - EXE' + B = 0   (Generalized Lyapunov Equation)
  32 ## @end group
  33 ## @end example
  34 ##
  35 ## @strong{Algorithm}@*
  36 ## Uses SLICOT SB03MD, SB04QD and SG03AD by courtesy of
  37 ## @uref{http://www.slicot.org, NICONET e.V.}
  38 ##
  39 ## @seealso{dlyapchol, lyap, lyapchol}
  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] = dlyap (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 ("dlyap: a, b must be real and square");
  55         error ("dlyap: %s, %s must be real and square", \
  56                 inputname (1), inputname (2));
  57       endif
  58
  59       if (rows (a) != rows (b))
  60         ## error ("dlyap: a, b must have the same number of rows");
  61         error ("dlyap: %s, %s must have the same number of rows", \
  62                 inputname (1), inputname (2));
  63       endif
  64
  65       [x, scale] = slsb03md (a, -b, true);     # AXA' - X = -B
  66
  67       ## x /= scale;                           # 0 < scale <= 1
  68
  69     case 3                                     # Sylvester equation
  70
  71       if (! is_real_square_matrix (a, b))
  72         ## error ("dlyap: a, b must be real and square");
  73         error ("dlyap: %s, %s must be real and square", \
  74                 inputname (1), inputname (2));
  75       endif
  76
  77       if (! is_real_matrix (c) || rows (c) != rows (a) || columns (c) != columns (b))
  78         ## error ("dlyap: c must be a real (%dx%d) matrix", rows (a), columns (b));
  79         error ("dlyap: %s must be a real (%dx%d) matrix", \
  80                 rows (a), columns (b), inputname (3));
  81       endif
  82
  83       x = slsb04qd (-a, b, c);                 # AXB' - X = -C
  84
  85     case 4                                     # generalized Lyapunov equation
  86
  87       if (! isempty (c))
  88         print_usage ();
  89       endif
  90
  91       if (! is_real_square_matrix (a, b, e))
  92         ## error ("dlyap: a, b, e must be real and square");
  93         error ("dlyap: %s, %s, %s must be real and square", \
  94                 inputname (1), inputname (2), inputname (4));
  95       endif
  96
  97       if (rows (b) != rows (a) || rows (e) != rows (a))
  98         ## error ("dlyap: a, b, e must have the same number of rows");
  99         error ("dlyap: %s, %s, %s must have the same number of rows", \
 100                 inputname (1), inputname (2), inputname (4));
 101       endif
 102
 103       if (! issymmetric (b))
 104         ## error ("dlyap: b must be symmetric");
 105         error ("dlyap: %s must be symmetric", \
 106                 inputname (2));
 107       endif
 108
 109       [x, scale] = slsg03ad (a, e, -b, true);  # AXA' - EXE' = -B
 110
 111       ## x /= scale;                           # 0 < scale <= 1
 112
 113     otherwise
 114       print_usage ();
 115
 116   endswitch
 117
 118   if (scale < 1)
 119     warning ("dlyap: solution scaled by %g to prevent overflow", scale);
 120   endif
 121
 122 endfunction
 123
 124
 125 ## Lyapunov
 126 %!shared X, X_exp
 127 %! A = [3.0   1.0   1.0
 128 %!      1.0   3.0   0.0
 129 %!      0.0   0.0   3.0];
 130 %!
 131 %! B = [25.0  24.0  15.0
 132 %!      24.0  32.0   8.0
 133 %!      15.0   8.0  40.0];
 134 %!
 135 %! X = dlyap (A.', -B);
 136 %!
 137 %! X_exp = [2.0000   1.0000   1.0000
 138 %!          1.0000   3.0000   0.0000
 139 %!          1.0000   0.0000   4.0000];
 140 %!
 141 %!assert (X, X_exp, 1e-4);
 142
 143 ## Sylvester
 144 %!shared X, X_exp
 145 %! A = [1.0   2.0   3.0
 146 %!      6.0   7.0   8.0
 147 %!      9.0   2.0   3.0];
 148 %!
 149 %! B = [7.0   2.0   3.0
 150 %!      2.0   1.0   2.0
 151 %!      3.0   4.0   1.0];
 152 %!
 153 %! C = [271.0   135.0   147.0
 154 %!      923.0   494.0   482.0
 155 %!      578.0   383.0   287.0];
 156 %!
 157 %! X = dlyap (-A, B, C);
 158 %!
 159 %! X_exp = [2.0000   3.0000   6.0000
 160 %!          4.0000   7.0000   1.0000
 161 %!          5.0000   3.0000   2.0000];
 162 %!
 163 %!assert (X, X_exp, 1e-4);
 164
 165 ## Generalized Lyapunov
 166 ## TODO: add a test