— 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

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.

— 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
     

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.

— 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
     

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

— Function File: dgram (a, b)

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
          

— Function File: dlyap (a, b)

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

— 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
     

or

              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.