diff scripts/control/base/lqr.m @ 3431:99ab64f4a09d

[project @ 2000-01-14 03:53:03 by jwe]
author jwe
date Fri, 14 Jan 2000 04:12:41 +0000
parents
children 22bd65326ec1
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/control/base/lqr.m	Fri Jan 14 04:12:41 2000 +0000
@@ -0,0 +1,174 @@
+## Copyright (C) 1993, 1994, 1995 Auburn University.  All rights reserved.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by the
+## Free Software Foundation; either version 2, or (at your option) any
+## later version.
+##
+## Octave is distributed in the hope that it will be useful, but WITHOUT
+## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+## FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+## for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{k}, @var{p}, @var{e}] =} lqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{z})
+## construct the linear quadratic regulator for the continuous time system
+## @iftex
+## @tex
+## $$
+##  {dx\over dt} = A x + B u
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+##
+## @example
+## dx
+## -- = A x + B u
+## dt
+## @end example
+##
+## @end ifinfo
+## to minimize the cost functional
+## @iftex
+## @tex
+## $$
+##  J = \int_0^\infty x^T Q x + u^T R u
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+##
+## @example
+##       infinity
+##       /
+##   J = |  x' Q x + u' R u
+##      /
+##     t=0
+## @end example
+## @end ifinfo
+##
+## @noindent
+## @var{z} omitted or
+## @iftex
+## @tex
+## $$
+##  J = \int_0^\infty x^T Q x + u^T R u + 2 x^T Z u
+## $$
+## @end tex
+## @end iftex
+## @ifinfo
+##
+## @example
+##       infinity
+##       /
+##   J = |  x' Q x + u' R u + 2 x' Z u
+##      /
+##     t=0
+## @end example
+##
+## @end ifinfo
+## @var{z} included.
+##
+## The following values are returned:
+##
+## @table @var
+## @item k
+## The state feedback gain,
+## @iftex
+## @tex
+## $(A - B K)$
+## @end tex
+## @end iftex
+## @ifinfo
+## (@var{a} - @var{b}@var{k})
+## @end ifinfo
+## is stable and minimizes the cost functional
+##
+## @item p
+## The stabilizing solution of appropriate algebraic Riccati equation.
+##
+## @item e
+## The vector of the closed loop poles of
+## @iftex
+## @tex
+## $(A - B K)$.
+## @end tex
+## @end iftex
+## @ifinfo
+## (@var{a} - @var{b}@var{k}).
+## @end ifinfo
+## @end table
+##
+## @strong{Reference}
+## Anderson and Moore, OPTIMAL CONTROL: LINEAR QUADRATIC METHODS,
+## Prentice-Hall, 1990, pp. 56-58
+## @end deftypefn
+
+## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu>
+## Created: August 1993.
+
+function [k, p, e] = lqr (a, b, q, r, s)
+
+  ## disp("lqr: entry");
+
+  if ((nargin != 4) && (nargin != 5))
+    error ("lqr: invalid number of arguments");
+  endif
+
+  ## Check a.
+  if ((n = is_square (a)) == 0)
+    error ("lqr: requires 1st parameter(a) to be square");
+  endif
+
+  ## Check b.
+  [n1, m] = size (b);
+  if (n1 != n)
+    error ("lqr: a,b not conformal");
+  endif
+
+  ## Check q.
+  if ( ((n1 = is_square (q)) == 0) || (n1 != n))
+    error ("lqr: q must be square and conformal with a");
+  endif
+
+  ## Check r.
+  if ( ((m1 = is_square(r)) == 0) || (m1 != m))
+    error ("lqr: r must be square and conformal with column dimension of b");
+  endif
+
+  ## Check if n is there.
+  if (nargin == 5)
+    [n1, m1] = size (s);
+    if ( (n1 != n) || (m1 != m))
+      error ("lqr: z must be identically dimensioned with b");
+    endif
+
+    ## Incorporate cross term into a and q.
+    ao = a - (b/r)*s';
+    qo = q - (s/r)*s';
+  else
+    s = zeros (n, m);
+    ao = a;
+    qo = q;
+  endif
+
+  ## Check that q, (r) are symmetric, positive (semi)definite
+
+  if (is_symmetric (q) && is_symmetric (r) ...
+      && all (eig (q) >= 0) && all (eig (r) > 0))
+    p = are (ao, (b/r)*b', qo);
+    k = r\(b'*p + s');
+    e = eig (a - b*k);
+  else
+    error ("lqr: q (r) must be symmetric positive (semi) definite");
+  endif
+
+  ## disp("lqr: exit");
+endfunction