@c Copyright (C) 1996, 1997 John W. Eaton
@c This is part of the Octave manual.
@c For copying conditions, see the file gpl.texi.

@node Control Theory, Signal Processing, Polynomial Manipulations, Top
@chapter Control Theory

Most of the functions described in this chapter were contributed by
A. Scottedward Hodel @email{A.S.Hodel@@eng.auburn.edu} and R. Bruce Tenison
@email{Bruce.Tenison@@eng.auburn.edu}.  They have also written a larger
collection of functions for solving linear control problems.  It is
currently being updated for Octave version 2, with snapshots of the
sources available from @url{ftp://ftp.eng.auburn.edu/pub/hodel}.

@deftypefn {Function File} {[@var{n}, @var{m}, @var{p}] =} abcddim (@var{a}, @var{b}, @var{c}, @var{d})
Check for compatibility of the dimensions of the matrices defining
the linear system
@iftex
@tex
$[A, B, C, D]$ corresponding to
$$
\eqalign{
 {dx\over dt} &= A x + B u\cr
            y &= C x + D u}
$$
@end tex
@end iftex
@ifinfo
[A, B, C, D] corresponding to

@example
dx/dt = a x + b u
y = c x + d u
@end example

@end ifinfo
or a similar discrete-time system.

If the matrices are compatibly dimensioned, then @code{abcddim} returns

@table @var
@item n
The number of system states.

@item m
The number of system inputs.

@item p
The number of system outputs.
@end table

Otherwise @code{abcddim} returns @var{n} = @var{m} = @var{p} = @minus{}1.
@end deftypefn

@deftypefn {Function File} {} are (@var{a}, @var{b}, @var{c}, @var{opt})

Return the solution, @var{x}, of the algebraic Riccati equation
@iftex
@tex
$$
A^TX + XA - XBX + C = 0
$$
@end tex
@end iftex
@ifinfo

@example
a' * x + x * a - x * b * x + c = 0
@end example
@end ifinfo

@noindent
for identically dimensioned square matrices @var{a}, @var{b}, and
@var{c}.  If @var{b} is not square, @code{are} attempts to use
@code{@var{b}*@var{b}'} instead.  If @var{c} is not square, @code{are}
attempts to use @code{@var{c}'*@var{c}}) instead.

To form the solution, Laub's Schur method (IEEE Transactions on
Automatic Control, 1979) is applied to the appropriate Hamiltonian
matrix.

The optional argument @var{opt} is passed to the eigenvalue balancing
routine.  If it is omitted, a value of @code{"B"} is assumed.
@end deftypefn

@deftypefn {Function File} {} c2d (@var{a}, @var{b}, @var{t})
Convert the continuous time system described by:
@iftex
@tex
$$
 {dx\over dt} = A x + B u
$$
@end tex
@end iftex
@ifinfo

@example
dx/dt = a x + b u
@end example

@end ifinfo
into a discrete time equivalent model
@iftex
@tex
$$
 x_{k+1} = A_d x_k + B_d u_k
$$
@end tex
@end iftex
@ifinfo

@example
x[k+1] = Ad x[k] + Bd u[k]
@end example
@end ifinfo

@noindent
via the matrix exponential assuming a zero-order hold on the input and
sample time @var{t}.
@end deftypefn

@deftypefn {Function File} {} dare (@var{a}, @var{b}, @var{c}, @var{r}, @var{opt})

Return the solution, @var{x} of the discrete-time algebraic Riccati
equation
@iftex
@tex
$$
A^TXA - X + A^TXB (R + B^TXB)^{-1} B^TXA + C = 0
$$
@end tex
@end iftex
@ifinfo

@example
a' x a - x + a' x b (r + b' x b)^(-1) b' x a + c = 0
@end example
@end ifinfo

@noindent
for matrices with dimensions:

@table @var
@item a
@var{n} by @var{n}.

@item b
@var{n} by @var{m}.

@item c
@var{n} by @var{n}, symmetric positive semidefinite.

@item r
@var{m} by @var{m}, symmetric positive definite (invertible).
@end table

If @var{c} is not square, then the function attempts to use
@code{@var{c}'*@var{c}} instead.

To form the solution, Laub's Schur method (IEEE Transactions on
Automatic Control, 1979) is applied to the appropriate symplectic
matrix.

See also Ran and Rodman, @cite{Stable Hermitian Solutions of Discrete
Algebraic Riccati Equations}, Mathematics of Control, Signals and 
Systems, Volume 5, Number 2 (1992).

The optional argument @var{opt} is passed to the eigenvalue balancing
routine.  If it is omitted, a value of @code{"B"} is assumed.
@end deftypefn

@deftypefn {Function File} {} dgram (@var{a}, @var{b})
Return the discrete controllability or observability gramian for the
discrete time system described by
@iftex
@tex
$$
 x_{k+1} = A x_k + B u_k
$$
$$
 y_k = C x_k + D u_k
$$
@end tex
@end iftex
@ifinfo

@example
x[k+1] = A x[k] + B u[k]
  y[k] = C x[k] + D u[k]
@end example
@end ifinfo

For example, @code{dgram (@var{a}, @var{b})} returns the discrete
controllability gramian and @code{dgram (@var{a}', @var{c}')} returns
the observability gramian.
@end deftypefn

@deftypefn {Function File} {[@var{l}, @var{m}, @var{p}, @var{e}] =} dlqe (@var{a}, @var{g}, @var{c}, @var{sigw}, @var{sigv}, @var{z})
Construct the linear quadratic estimator (Kalman filter) for the
discrete time system
@iftex
@tex
$$
 x_{k+1} = A x_k + B u_k + G w_k
$$
$$
 y_k = C x_k + D u_k + w_k
$$
@end tex
@end iftex
@ifinfo

@example
x[k+1] = A x[k] + B u[k] + G w[k]
  y[k] = C x[k] + D u[k] + w[k]
@end example

@end ifinfo
where @var{w}, @var{v} are zero-mean gaussian noise processes with
respective intensities @code{@var{sigw} = cov (@var{w}, @var{w})} and
@code{@var{sigv} = cov (@var{v}, @var{v})}.

If specified, @var{z} is @code{cov (@var{w}, @var{v})}.  Otherwise
@code{cov (@var{w}, @var{v}) = 0}.

The observer structure is
@iftex
@tex
$$
 z_{k+1} = A z_k + B u_k + k(y_k - C z_k - D u_k)
$$
@end tex
@end iftex
@ifinfo

@example
z[k+1] = A z[k] + B u[k] + k(y[k] - C z[k] - D u[k])
@end example
@end ifinfo

@noindent
The following values are returned:

@table @var
@item l
The observer gain, 
@iftex
@tex
$(A - ALC)$.
@end tex
@end iftex
@ifinfo
(@var{a} - @var{a}@var{l}@var{c}).
@end ifinfo
is stable.

@item m
The Riccati equation solution.

@item p
The estimate error covariance after the measurement update.

@item e
The closed loop poles of
@iftex
@tex
$(A - ALC)$.
@end tex
@end iftex
@ifinfo
(@var{a} - @var{a}@var{l}@var{c}).
@end ifinfo
@end table
@end deftypefn

@deftypefn {Function File} {[@var{k}, @var{p}, @var{e}] =} dlqr (@var{a}, @var{b}, @var{q}, @var{r}, @var{z})
Construct the linear quadratic regulator for the discrete time system
@iftex
@tex
$$
 x_{k+1} = A x_k + B u_k
$$
@end tex
@end iftex
@ifinfo

@example
x[k+1] = A x[k] + B u[k]
@end example

@end ifinfo
to minimize the cost functional
@iftex
@tex
$$
 J = \sum x^T Q x + u^T R u
$$
@end tex
@end iftex
@ifinfo

@example
J = Sum (x' Q x + u' R u)
@end example
@end ifinfo

@noindent
@var{z} omitted or
@iftex
@tex
$$
 J = \sum x^T Q x + u^T R u + 2 x^T Z u
$$
@end tex
@end iftex
@ifinfo

@example
J = Sum (x' Q x + u' R u + 2 x' Z u)
@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.

@item p
The solution of algebraic Riccati equation.

@item e
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
@end deftypefn

@deftypefn {Function File} {} dlyap (@var{a}, @var{b})
Solve the discrete-time Lyapunov equation
@iftex
@tex
$AXA^T - X + B = 0$
@end tex
@end iftex
@ifinfo
@code{a x a' - x + b = 0}
@end ifinfo
for square matrices @var{a}, @var{b}.  If @var{b} is not square, then the
function attempts to solve either
@iftex
@tex
$AXA^T - X + B B^T = 0$ or $A^TXA - X + B^TB = 0$,
@end tex
@ifinfo
@code{a x a' - x + b b' = 0} or @code{a' x a - x + b' b = 0},
@end ifinfo
whichever is appropriate.

Uses Schur decomposition method as in Kitagawa 
@iftex
@tex
@cite{An Algorithm for Solving the Matrix Equation $X = F X F^\prime + S$},
@end tex
@end iftex
@ifinfo
@cite{An Algorithm for Solving the Matrix Equation @var{X} =
@var{F}@var{X}@var{F}' + @var{S}},
@end ifinfo
International Journal of Control, Volume 25, Number 5, pages 745--753
(1977); column-by-column solution method as suggested in
Hammerling, @cite{Numerical Solution of the Stable, Non-Negative
Definite Lyapunov Equation}, IMA Journal of Numerical Analysis, Volume
2, pages 303--323 (1982).
@end deftypefn

@deftypefn {Function File} {} is_controllable (@var{a}, @var{b}, @var{tol})
Return 1 if the pair (@var{a}, @var{b}) is controllable.  Otherwise,
return 0.

The optional argument @var{tol} is a roundoff parameter.  If it is
omitted, a value of @code{2*eps} is used.

Currently, @code{is_controllable} just constructs the controllability
matrix and checks rank.
@end deftypefn

@deftypefn {Function File} {} is_observable (@var{a}, @var{c}, @var{tol})

Return 1 if the pair (@var{a}, @var{c}) is observable.
Otherwise, return 0.

The optional argument @var{tol} is a roundoff parameter.  If it is
omitted, a value of @code{2*eps} is used.
@end deftypefn

@deftypefn {Function File} {[@var{k}, @var{p}, @var{e}] =} lqe (@var{a}, @var{g}, @var{c}, @var{sigw}, @var{sigv}, @var{z})
Construct the linear quadratic estimator (Kalman filter) for the
continuous time system
@iftex
@tex
$$
 {dx\over dt} = A x + B u
$$
$$
 y = C x + D u
$$
@end tex
@end iftex
@ifinfo

@example
dx
-- = a x + b u
dt

y = c x + d u
@end example

@end ifinfo
where @var{w} and @var{v} are zero-mean gaussian noise processes with
respective intensities

@example
sigw = cov (w, w)
sigv = cov (v, v)
@end example

The optional argument @var{z} is the cross-covariance
@code{cov (@var{w}, @var{v})}.  If it is omitted,
@code{cov (@var{w}, @var{v}) = 0} is assumed.

Observer structure is @code{dz/dt = A z + B u + k (y - C z - D u)}

The following values are returned:

@table @var
@item k
The observer gain,
@iftex
@tex
$(A - K C)$
@end tex
@end iftex
@ifinfo
(@var{a} - @var{k}@var{c})
@end ifinfo
is stable.

@item p
The solution of algebraic Riccati equation.

@item e
The vector of closed loop poles of
@iftex
@tex
$(A - K C)$.
@end tex
@end iftex
@ifinfo
(@var{a} - @var{k}@var{c}).
@end ifinfo
@end table
@end deftypefn

@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.

@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
@end deftypefn

@deftypefn {Function File} {} lyap (@var{a}, @var{b}, @var{c})
Solve the Lyapunov (or Sylvester) equation via the Bartels-Stewart
algorithm (Communications of the ACM, 1972).

If @var{a}, @var{b}, and @var{c} are specified, then @code{lyap} returns
the solution of the  Sylvester equation
@iftex
@tex
$$
 A X + X B + C = 0
$$
@end tex
@end iftex
@ifinfo

@example
a x + x b + c = 0
@end example
@end ifinfo

If only @code{(a, b)} are specified, then @code{lyap} returns the
solution of the Lyapunov equation
@iftex
@tex
$$
 A^T X + X A + B = 0
$$
@end tex
@end iftex
@ifinfo

@example
a' x + x a + b = 0
@end example
@end ifinfo

If @var{b} is not square, then @code{lyap} returns the solution of either
@iftex
@tex
$$
 A^T X + X A + B^T B = 0
$$
@end tex
@end iftex
@ifinfo

@example
a' x + x a + b' b = 0
@end example
@end ifinfo

@noindent
or
@iftex
@tex
$$
 A X + X A^T + B B^T = 0
$$
@end tex
@end iftex
@ifinfo

@example
a x + x a' + b b' = 0
@end example
@end ifinfo

@noindent
whichever is appropriate.
@end deftypefn

@deftypefn {Function File} {} tzero (@var{a}, @var{b}, @var{c}, @var{d}, @var{opt})
Compute the transmission zeros of
@iftex
@tex
$[A, B, C, D]$.
@end tex
@end iftex
@ifinfo
[A, B, C, D].
@end ifinfo

The optional argument @var{opt} is passed to the eigenvalue balancing
routine.  If it is omitted, a value of @code{"B"} is assumed.
@end deftypefn