octave_packages/control-2.3.52/place.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{f} =} place (@var{sys}, @var{p})
  20 ## @deftypefnx {Function File} {@var{f} =} place (@var{a}, @var{b}, @var{p})
  21 ## @deftypefnx {Function File} {[@var{f}, @var{info}] =} place (@var{sys}, @var{p}, @var{alpha})
  22 ## @deftypefnx {Function File} {[@var{f}, @var{info}] =} place (@var{a}, @var{b}, @var{p}, @var{alpha})
  23 ## Pole assignment for a given matrix pair (@var{A},@var{B}) such that @code{p = eig (A-B*F)}.
  24 ## If parameter @var{alpha} is specified, poles with real parts (continuous-time)
  25 ## or moduli (discrete-time) below @var{alpha} are left untouched.
  26 ##
  27 ## @strong{Inputs}
  28 ## @table @var
  29 ## @item sys
  30 ## LTI system.
  31 ## @item a
  32 ## State transition matrix (n-by-n) of a continuous-time system.
  33 ## @item b
  34 ## Input matrix (n-by-m) of a continuous-time system.
  35 ## @item p
  36 ## Desired eigenvalues of the closed-loop system state-matrix @var{A-B*F}.
  37 ## @code{length (p) <= rows (A)}.
  38 ## @item alpha
  39 ## Specifies the maximum admissible value, either for real
  40 ## parts or for moduli, of the eigenvalues of @var{A} which will
  41 ## not be modified by the eigenvalue assignment algorithm.
  42 ## @code{alpha >= 0} for discrete-time systems.
  43 ## @end table
  44 ##
  45 ## @strong{Outputs}
  46 ## @table @var
  47 ## @item f
  48 ## State feedback gain matrix.
  49 ## @item info
  50 ## Structure containing additional information.
  51 ## @item info.nfp
  52 ## The number of fixed poles, i.e. eigenvalues of @var{A} having
  53 ## real parts less than @var{alpha}, or moduli less than @var{alpha}.
  54 ## These eigenvalues are not modified by @command{place}.
  55 ## @item info.nap
  56 ## The number of assigned eigenvalues.  @code{nap = n-nfp-nup}.
  57 ## @item info.nup
  58 ## The number of uncontrollable eigenvalues detected by the
  59 ## eigenvalue assignment algorithm.
  60 ## @item info.z
  61 ## The orthogonal matrix @var{z} reduces the closed-loop
  62 ## system state matrix @code{A + B*F} to upper real Schur form.
  63 ## Note the positive sign in @code{A + B*F}.
  64 ## @end table
  65 ##
  66 ## @strong{Note}
  67 ## @example
  68 ## Place is also suitable to design estimator gains:
  69 ## @group
  70 ## L = place (A.', C.', p).'
  71 ## L = place (sys.', p).'   # useful for discrete-time systems
  72 ## @end group
  73 ## @end example
  74 ##
  75 ## @strong{Algorithm}@*
  76 ## Uses SLICOT SB01BD by courtesy of
  77 ## @uref{http://www.slicot.org, NICONET e.V.}
  78 ## @end deftypefn
  79
  80 ## Special thanks to Peter Benner from TU Chemnitz for his advice.
  81 ## Author: Lukas Reichlin <lukas.reichlin@gmail.com>
  82 ## Created: December 2009
  83 ## Version: 0.5
  84
  85 function [f, info] = place (a, b, p = [], alpha = [], tol = [])
  86
  87   if (nargin < 2 || nargin > 5)
  88     print_usage ();
  89   endif
  90
  91   if (isa (a, "lti"))              # place (sys, p), place (sys, p, alpha), place (sys, p, alpha, tol)
  92     if (nargin > 4)                # nargin < 2 already tested
  93       print_usage ();
  94     endif
  95     tol = alpha;
  96     alpha = p;
  97     p = b;
  98     sys = a;
  99     [a, b] = ssdata (sys);         # descriptor matrice e should be regular
 100     discrete = ! isct (sys);       # treat tsam = -2 as continuous system
 101   else                             # place (a, b, p), place (a, b, p, alpha), place (a, b, p, alpha, tol)
 102     if (nargin < 3)                # nargin > 5 already tested
 103       print_usage ();
 104     endif
 105     if (! is_real_square_matrix (a) || ! is_real_matrix (b) || rows (a) != rows (b))
 106       error ("place: matrices a and b not conformal");
 107     endif
 108     discrete = 0;                  # assume continuous system
 109   endif
 110
 111   if (! isnumeric (p) || ! isvector (p) || isempty (p))  # p could be complex
 112     error ("place: p must be a vector");
 113   endif
 114
 115   p = sort (reshape (p, [], 1));   # complex conjugate pairs must appear together
 116   wr = real (p);
 117   wi = imag (p);
 118
 119   n = rows (a);                    # number of states
 120   np = length (p);                 # number of given eigenvalues
 121
 122   if (np > n)
 123     error ("place: at most %d eigenvalues can be assigned for the given matrix a (%dx%d)",
 124             n, n, n);
 125   endif
 126
 127   if (isempty (alpha))
 128     if (discrete)
 129       alpha = 0;
 130     else
 131       alpha = - norm (a, inf);
 132     endif
 133   endif
 134
 135   if (isempty (tol))
 136     tol = 0;
 137   endif
 138
 139   [f, nfp, nap, nup, z] = slsb01bd (a, b, wr, wi, discrete, alpha, tol);
 140   f = -f;                          # A + B*F --> A - B*F
 141
 142   info = struct ("nfp", nfp, "nap", nap, "nup", nup, "z", z);
 143
 144 endfunction
 145
 146
 147 ## Test from "legacy" control package 1.0.*
 148 %!shared A, B, C, P, Kexpected
 149 %! A = [0, 1; 3, 2];
 150 %! B = [0; 1];
 151 %! C = [2, 1];  # C is needed for ss; it doesn't matter what the value of C is
 152 %! P = [-1, -0.5];
 153 %! Kexpected = [3.5, 3.5];
 154 %!assert (place (ss (A, B, C), P), Kexpected, 2*eps);
 155 %!assert (place (A, B, P), Kexpected, 2*eps);
 156
 157 ## FIXME: Test from SLICOT example SB01BD fails with 4 eigenvalues in P
 158 %!shared F, F_exp, ev_ol, ev_cl
 159 %! A = [-6.8000   0.0000  -207.0000   0.0000
 160 %!       1.0000   0.0000     0.0000   0.0000
 161 %!      43.2000   0.0000     0.0000  -4.2000
 162 %!       0.0000   0.0000     1.0000   0.0000];
 163 %!
 164 %! B = [ 5.6400   0.0000
 165 %!       0.0000   0.0000
 166 %!       0.0000   1.1800
 167 %!       0.0000   0.0000];
 168 %!
 169 %! P = [-0.5000 + 0.1500i
 170 %!      -0.5000 - 0.1500i];
 171 #%!      -2.0000 + 0.0000i
 172 #%!      -0.4000 + 0.0000i];
 173 %!
 174 %! ALPHA = -0.4;
 175 %! TOL = 1e-8;
 176 %!
 177 %! F = place (A, B, P, ALPHA, TOL);
 178 %!
 179 %! F_exp = - [-0.0876  -4.2138   0.0837 -18.1412
 180 %!            -0.0233  18.2483  -0.4259  -4.8120];
 181 %!
 182 %! ev_ol = sort (eig (A));
 183 %! ev_cl = sort (eig (A - B*F));
 184 %!
 185 %!assert (F, F_exp, 1e-4);
 186 %!assert (ev_ol(3:4), ev_cl(3:4), 1e-4);