Periodic solver utilities

 pssylvdckr(A, B, C; adj = true, isgn = 1) -> X

For the periodic matrices A, B and C compute the periodic matrix X, which satisfies the periodic discrete-time Sylvester matrix equation of continuous-time flavour

  isgn*A(i)*X(i+1) + X(i)*B(i) = C(i), if rev = true,

or

  A(i)*X(i) + isgn*X(i+1)*B(i) =  C(i), if rev = false,

where abs(isgn) = 1. The periodic matrices A, B and C have periods pa, pb and pc, respectively, and are stored as either 3-dimensional arrays or as vectors of matrices. The resulting periodic matrix X has period p = LCM(pa,pb,pc) and is stored correspondingly.

For the solvability of the above equations the following conditions must hold: if rev = false, the matrix products A(p)*A(p-1)...A(1) and B(p)*B(p-1)*...B(1) must not have eigenvalues α and β such that α +(-isgn)^p*β = 0, while if rev = true, the matrix products A(1)*A(2)...A(p) and B(1)*B(2)*...B(p) must not have eigenvalues α and β such that α +(-isgn)^p*β = 0.

The Kronecker product expansion of equations is employed and therefore this function is not recommended for large order matrices or large periods.

pssylvdc(A, B, C; rev = false, isgn = 1, fast = true) -> X

Solve the periodic discrete-time Sylvester equation of continuous-time flavour.

For the square n-th order periodic matrices A(i), i = 1, ..., pa, B(i), i = 1, ..., pb and C(i), i = 1, ..., pc of periods pa, pb and pc, respectively, the periodic solution X(i), i = 1, ..., p of period p = lcm(pa,pb,pc) of the periodic Sylvester equation is computed:

A(i)*X(i) + isgn*X(i+1)*B(i) = C(i), i = 1, ..., p     for `rev = false`; 

isgn*A(i)*X(i+1) + X(i)*B(i) = C(i), i = 1, ..., p     for `rev = true`.

The periodic matrices A, B and C are stored in the m×m×pa, n×n×pb and m×n×pc 3-dimensional arrays A, B and C, respectively, and X results as a m×n×p 3-dimensional array.

Alternatively, the periodic matrices A, B and C can be stored in the pa-, pb- and pc-dimensional vectors of matrices A, B and C, respectively, and X results as a p-dimensional vector of matrices.

The periodic discrete analog of the Bartels-Stewart method based on the periodic Schur form of the periodic matrices A and B is employed (see Appendix II of [1]). If fast = true, the QR factorization of bordered-almost-block-diagonal (BABD) matrix algorithm of [2] is employed to solve periodic Sylvester equations up to order 2. This option is more appropriate for large periods. If fast = false, the QR factorization of the cyclic Kronecker form for the periodic Sylvester operator is used to to solve periodic Sylvester equations up to order 2.

For the existence of a solution A and B must not have characteristic multipliers α and β such that α +isgn*β = 0.

Reference:

[1] A. Varga. Robust and minimum norm pole assignment with periodic state feedback. IEEE Trans. on Automatic Control, vol. 45, pp. 1017-1022, 2000.

[2] R. Granat, B. Kågström, and D. Kressner, Computing periodic deflating subspaces associated with a specified set of eigenvalues. BIT Numerical Mathematics vol. 47, pp. 763–791, 2007.

pslyapd(A, C; adj = true, stability_check = false) -> X

Solve the periodic discrete-time Lyapunov equation.

For the square n-th order periodic matrices A(i), i = 1, ..., pa and C(i), i = 1, ..., pc of periods pa and pc, respectively, the periodic solution X(i), i = 1, ..., p of period p = lcm(pa,pc) of the periodic Lyapunov equation is computed:

A(i)'*X(i+1)*A(i) + C(i) = X(i), i = 1, ..., p     for `adj = true`; 

A(i)*X(i)*A(i)' + C(i) = X(i+1), i = 1, ..., p     for `adj = false`.

The periodic matrices A and C are stored in the n×n×pa and n×n×pc 3-dimensional arrays A and C, respectively, and X results as a n×n×p 3-dimensional array.

Alternatively, the periodic matrices A and C can be stored in the pa- and pc-dimensional vectors of matrices A and C, respectively, and X results as a p-dimensional vector of matrices.

If stability_check = true, the stability of characteristic multipliers of A is checked and an error is issued if any characteristic multiplier has modulus equal to or larger than one.

The periodic discrete analog of the Bartels-Stewart method based on the periodic Schur form of the periodic matrix A is employed [1].

Reference:

[1] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol. 67, pp. 69-87, 1997.

pslyapd2(A, C, E; stability_check = false) -> (X, Y)

Solve a pair of periodic discrete-time Lyapunov equations.

For the square n-th order periodic matrices A(i), i = 1, ..., pa, C(i), i = 1, ..., pc, and E(i), i = 1, ..., pe of periods pa, pc and pe, respectively, the periodic solutions X(i), i = 1, ..., p and Y(i), i = 1, ..., p of period p = lcm(pa,pc,pe) of the periodic Lyapunov equations are computed:

A(i)*X(i)*A(i)' + C(i) = X(i+1), i = 1, ..., p ,  
 
A(i)'*Y(i+1)*A(i) + E(i) = Y(i), i = 1, ..., p .

The periodic matrices A, C and E are stored in the n×n×pa, n×n×pc and n×n×pe 3-dimensional arrays A, C and E, respectively, and X and Y result as n×n×p 3-dimensional arrays.

Alternatively, the periodic matrices A, C and E can be stored in the pa-, pc- and pe-dimensional vectors of matrices A, C and E, respectively, and X and Y result as p-dimensional vectors of matrices.

If stability_check = true, the stability of characteristic multipliers of A is checked and an error is issued if any characteristic multiplier has modulus equal to or larger than one.

The periodic discrete analog of the Bartels-Stewart method based on the periodic Schur form of the periodic matrix A is employed [1].

Reference:

[1] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol. 67, pp. 69-87, 1997.

 pslyapdkr(A, C; adj = true) -> X

Solve the periodic discrete-time Lyapunov matrix equation

  A'σXA + C = X, if adj = true,

or

  A*X*A' + C =  σX, if adj = false,

where σ is the forward shift operator σX(i) = X(i+1). The periodic matrix A must not have characteristic multipliers on the unit circle. The periodic matrices A and C are either stored as 3-dimensional arrays or as as vectors of matrices.

The Kronecker product expansion of equations is employed and therefore this function is not recommended for large order matrices or large periods.

 psplyapd(A, C; adj = true, rtol = ϵ^(3/4)) -> U

Compute the upper triangular factor U of the solution X = U'U of the periodic discrete-time Lyapunov matrix equation

  A'σXA + C'C = X, if adj = true,

or of the solution X = UU' of the periodic discrete-time Lyapunov matrix equation

  AXA' + CC' =  σX, if adj = false,

where σ is the forward shift operator σX(i) = X(i+1). The periodic matrix A must be stable, i.e., have all characteristic multipliers with moduli less than one.

The periodic matrices A and C are either stored as 3-dimensional arrays or as as vectors of matrices.

The iterative method (Algorithm 5) of [1] and its dual version are employed.

Reference:

[1] A. Varga, "Periodic Lyapunov equations: some applications and new algorithms", Int. J. Control, vol. 67, pp. 69-87, 1997.