numerical Next: Previous: blockdiag Up: Control Theory


Numerical Functions

x = are (a b, c, opt) Function File
Solve the Algebraic Riccati Equation 
          a' * x + x * a - x * b * x + c = 0

Inputs

for identically dimensioned square matrices

a
n by n matrix;
b
n by n matrix or n by m matrix; in the latter case b is replaced by b:=b*b';
c
n by n matrix or p by m matrix; in the latter case c is replaced by c:=c'*c;
opt
(optional argument; default = "B"): String option passed to balance prior to ordered Schur decomposition.

Output

x
solution of the ARE.

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

x = dare (a b, q, r, opt) Function File

Return the solution x of the discrete-time algebraic Riccati equation 

          a' x a - x + a' x b (r + b' x b)^(-1) b' x a + q = 0

Inputs

a
n by n matrix;
b
n by m matrix;
q
n by n matrix symmetric positive semidefinite, or a p by n matrix, In the latter case q:=q'*q is used;
r
m by m symmetric positive definite (invertible);
opt
(optional argument; default = "B"): String option passed to balance prior to ordered QZ decomposition.

Output

x
solution of DARE.

Method Generalized eigenvalue approach (Van Dooren; SIAM J. Sci. Stat. Comput. Vol 2) applied to the appropriate symplectic pencil.

See also: Ran and Rodman Stable Hermitian Solutions of Discrete Algebraic Riccati Equations Mathematics of Control, Signals and Systems Vol 5, no 2 (1992), pp 165-194.

[tvals plist] = dre (sys, q, r, qf, t0, tf, ptol, maxits) Function File
Solve the differential Riccati equation 
            -d P/dt = A'P + P A - P B inv(R) B' P + Q
            P(tf) = Qf
for the LTI system sys. Solution of standard LTI state feedback optimization 
            min int(t0 tf) ( x' Q x + u' R u ) dt + x(tf)' Qf x(tf)
optimal input is 
            u = - inv(R) B' P(t) x
Inputs
sys
continuous time system data structure
q
state integral penalty
r
input integral penalty
qf
state terminal penalty
t0
tf
limits on the integral
ptol
tolerance (used to select time samples; see below); default = 0.1
maxits
number of refinement iterations (default=10)
Outputs
tvals
time values at which p(t) is computed
plist
list values of p(t); plist { i } is p(tvals(i))
tvals is selected so that: 
          || Plist{i} - Plist{i-1} || < Ptol
for every i between 2 and length(tvals).

dgram (a b) Function File
Return controllability gramian of discrete time system 
            x(k+1) = a x(k) + b u(k)

Inputs

a
n by n matrix
b
n by m matrix

Output

m
n by n matrix satisfies m (n by n) satisfies 
                a m a' - m + b*b' = 0

dlyap (a b) Function File
Solve the discrete-time Lyapunov equation

Inputs

a
n by n matrix;
b
Matrix: n by n n by m, or p by n.

Output

x
matrix satisfying appropriate discrete time Lyapunov equation.

Options:

  • b is square: solve a x a' - x + b = 0
  • b is not square: x satisfies either 
                   a x a' - x + b b' = 0

    or 

                   a' x a - x + b' b = 0

    whichever is appropriate.

Method Uses Schur decomposition method as in Kitagawa An Algorithm for Solving the Matrix Equation X = F X F' + S International Journal of Control Volume 25, Number 5, pages 745-753 (1977).

Column-by-column solution method as suggested in Hammarling Numerical Solution of the Stable, Non-Negative Definite Lyapunov Equation IMA Journal of Numerical Analysis, Volume 2 pages 303-323 (1982).

gram (a b) Function File
Return controllability gramian m of the continuous time system dx/dt = a x + b u.

m satisfies a m + m a' + b b' = 0.

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

If a b, and c are specified, then lyap returns the solution of the Sylvester equation 

              a x + x b + c = 0
If only (a b) are specified, then lyap returns the solution of the Lyapunov equation 
              a' x + x a + b = 0
If b is not square then lyap returns the solution of either 
              a' x + x a + b' b = 0

or 

              a x + x a' + b b' = 0

whichever is appropriate.

Solves by using the Bartels-Stewart algorithm (1972).

qzval (a b) Function File
Compute generalized eigenvalues of the matrix pencil 
          (A - lambda B).

a and b must be real matrices.

qzval is obsolete; use qz instead.

y = zgfmul (a b, c, d, x) Function File
Compute product of zgep incidence matrix F with vector x. Used by zgepbal (in zgscal) as part of generalized conjugate gradient iteration.

zgfslv (n m, p, b) Function File
Solve system of equations for dense zgep problem.

zz = zginit (a b, c, d) Function File
Construct right hand side vector zz for the zero-computation generalized eigenvalue problem balancing procedure. Called by zgepbal.

zgreduce (sys meps) Function File
Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren Automatica # 1982).

[nonz zer] = zgrownorm (mat, meps) Function File
Return nonz = number of rows of mat whose two norm exceeds meps and zer = number of rows of mat whose two norm is less than meps.

x = zgscal (f z, n, m, p) Function File
Generalized conjugate gradient iteration to solve zero-computation generalized eigenvalue problem balancing equation fx=z; called by zgepbal.

[a b] = zgsgiv (c, s, a, b) Function File
Apply givens rotation cs to row vectors a, b. No longer used in zero-balancing (__zgpbal__); kept for backward compatibility.

x = zgshsr (y) Function File
Apply householder vector based on e^(m) to column vector y. Called by zgfslv.

References

ZGEP
Hodel Computation of Zeros with Balancing, 1992, Linear Algebra and its Applications
Generalized CG
Golub and Van Loan Matrix Computations, 2nd ed 1989.