Next: , Previous: blockdiag, Up: Control Theory

29.5 Numerical Functions

— Function File: x = are (a, b, c, opt)

Solve the Algebraic Riccati Equation

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

Inputs for identically dimensioned square matrices

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


solution of the ARE.

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

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

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


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


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.

— Function File: [tvals, plist] = dre (sys, q, r, qf, t0, tf, ptol, maxits)

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


continuous time system data structure
state integral penalty
input integral penalty
state terminal penalty
limits on the integral
tolerance (used to select time samples; see below); default = 0.1
number of refinement iterations (default=10)
time values at which p(t) is computed
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).

— Function File: dgram (a, b)

Return controllability gramian of discrete time system

            x(k+1) = a x(k) + b u(k)


n by n matrix
n by m matrix


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

— Function File: dlyap (a, b)

Solve the discrete-time Lyapunov equation


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


matrix satisfying appropriate discrete time Lyapunov equation.


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

— Function File: gram (a, b)

Return controllability gramian m of the continuous time system dx/dt = a x + b u.

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

— Function File: lyap (a, b, c)
— Function File: lyap (a, b)

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


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

whichever is appropriate.

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

— Function File: qzval (a, b)

Compute generalized eigenvalues of the matrix pencil

          (A - lambda B).

a and b must be real matrices.

qzval is obsolete; use qz instead.

— Function File: y = zgfmul (a, b, c, d, x)

Compute product of zgep incidence matrix F with vector x. Used by zgepbal (in zgscal) as part of generalized conjugate gradient iteration.

— Function File: zgfslv (n, m, p, b)

Solve system of equations for dense zgep problem.

— Function File: zz = zginit (a, b, c, d)

Construct right hand side vector zz for the zero-computation generalized eigenvalue problem balancing procedure. Called by zgepbal.

— Function File: zgreduce (sys, meps)

Implementation of procedure REDUCE in (Emami-Naeini and Van Dooren, Automatica, # 1982).

— Function File: [nonz, zer] = zgrownorm (mat, meps)

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.

— Function File: x = zgscal (f, z, n, m, p)

Generalized conjugate gradient iteration to solve zero-computation generalized eigenvalue problem balancing equation fx=z; called by zgepbal.

— Function File: [a, b] = zgsgiv (c, s, a, b)

Apply givens rotation c,s to row vectors a, b. No longer used in zero-balancing (__zgpbal__); kept for backward compatibility.

— Function File: x = zgshsr (y)

Apply householder vector based on e^(m) to column vector y. Called by zgfslv.


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