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

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.

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.