Periodic Sylvester equation solvers

pcsylv(A, C; K = 10, N = 1, adj = false, solver, reltol, abstol, intpol = true) -> X

Solve the periodic Sylvester differential equation

.
X(t) = A(t)X(t) + X(t)B(t) + C(t) , if adj = false,

or

 .
-X(t) = A(t)X(t) + X(t)B(t) + C(t) , if adj = true.

The periodic matrices A, B and C must have the same type, compatible dimensions and commensurate periods. Additionally, A and B must be square. The resulting periodic solution X has the type PeriodicFunctionMatrix and X(t) can be used to evaluate the value of X at time t. X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

The multiple-shooting method of [1] is employed to convert the (continuous-time) periodic differential Sylvester equation into a discrete-time periodic Sylvester equation satisfied by a multiple point periodic generator of the solution. The keyword argument K specifies the number of grid points to be used for the discretization of the continuous-time problem (default: K = 10). If A, B and C are of types PeriodicTimeSeriesMatrix or PeriodicSwitchingMatrix, then K specifies the number of grid points used between two consecutive switching time values (default: K = 1). The multiple point periodic generator is computed by solving the appropriate discrete-time periodic Sylvester equation using the periodic Schur method of [2]. The obtained periodic generator is finally converted into a periodic function matrix which determines for a given t the function value X(t) by interpolating the solution of the appropriate differential equations if intpol = true (default) or by integrating the appropriate ODE from the nearest grid point value if intpol = false. For the interpolation-based evaluation the integer keyword argument N can be used to split the integration domain (i.e., one period) into N subdomains to perform the interpolations separately in each domain. The default value of N is N = 1.

The ODE solver to be employed to convert the continuous-time problem into a discrete-time problem can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4) and absolute accuracy abstol (default: abstol = 1.e-7). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

Parallel computation of the matrices of the discrete-time problem can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

References

[1] A. Varga. On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. IEEE CDC/ECC, Seville, 2005.

[2] R. Byers and N. Rhee, “Cyclic Schur and Hessenberg-Schur numerical methods for solving periodic Sylvester and Sylvester equations,” Dept. of Mathematics, Univ. of Missouri at Kansas City, Tech. Rep., June 1995.

prcsylv(A, B, C; K = 10, N = 1, solver, reltol, abstol, intpol) -> X

Solve the periodic reverse-time Sylvester differential equation

 .
-X(t) = A(t)'X(t) + X(t)A(t) + C(t).

The periodic matrices A, B and C must have the same type, compatible dimensions and commensurate periods. Additionally, A and B must be square. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

This function is merely an interface to pcsylv (see this function for the description of keyword parameters).

pfcsylv(A, B, C; K = 10, N = 1, solver, reltol, abstol, intpol) -> X

Solve the periodic forward-time Sylvester differential equation

.
X(t) = A(t)X(t) + X(t)B(t) + C(t) .

The periodic matrices A, B and C must have the same type, compatible dimensions and commensurate periods, Additionally, A and B must be square. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

This function is merely an interface to pcsylv (see this function for the description of keyword parameters).

pgsylv(A, B, C[, K = 1]; adj = false, solver, reltol, abstol) -> X

Compute periodic generators for the periodic Sylvester differential equation

.
X(t) = A(t)X(t) + X(t)B(t) + C(t) , if adj = false,

or

 .
-X(t) = A(t)X(t) + X(t)B(t) + C(t) , if adj = true.

The periodic matrices A, B and C must have the same type, compatible dimensions and commensurate periods. K specifies the number of grid points to be used for the discretization of the continuous-time problem (default: K = 1). If A, B and C are of types PeriodicTimeSeriesMatrix or PeriodicSwitchingMatrix, then K specifies the number of grid points used between two consecutive switching time values (default: K = 1).

If A, B and C have the types PeriodicFunctionMatrix, HarmonicArray, FourierFunctionMatrix or PeriodicTimeSeriesMatrix, then the resulting X is a collection of periodic generator matrices determined as a periodic time-series matrix with N components, where N = 1 if A, B and C are constant matrices and N = K otherwise. If A, B and C have the type PeriodicSwitchingMatrix, then X is a collection of periodic generator matrices determined as a periodic switching matrix, whose switching times are the unique switching times contained in the union of the switching times of A, B and C. If K > 1, a refined grid of K equidistant values is used for each two consecutive switching times in the union.

The period of X is set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly. Any component matrix of X is a valid initial value to be used to generate the solution over a full period by integrating the appropriate differential equation. The multiple-shooting method of [1] is employed, first, to convert the continuous-time periodic Sylvester differential equation into a discrete-time periodic Sylvester equation satisfied by the generator solution in the grid points and then to compute the solution by solving an appropriate discrete-time periodic Sylvester equation using the periodic Schur method of [2].

The ODE solver to be employed to convert the continuous-time problem into a discrete-time problem can be specified using the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4) and absolute accuracy abstol (default: abstol = 1.e-7). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

Parallel computation of the matrices of the discrete-time problem can be alternatively performed by starting Julia with several execution threads. The number of execution threads is controlled either by using the -t/--threads command line argument or by using the JULIA_NUM_THREADS environment variable.

References

[1] A. Varga. On solving periodic differential matrix equations with applications to periodic system norms computation. Proc. IEEE CDC/ECC, Seville, 2005.

[2] R. Byers and N. Rhee, “Cyclic Schur and Hessenberg-Schur numerical methods for solving periodic Sylvester and Sylvester equations,” Dept. of Mathematics, Univ. of Missouri at Kansas City, Tech. Rep., June 1995.

  tvcsylv_eval(t, W, A, C; adj = false, solver, reltol, abstol) -> Xval

Compute the time value Xval := X(t) of the solution of the periodic Sylvester differential equation

   .
   X(t) = A(t)X(t) + X(t)B(t) + C(t) ,  X(t0) = W(t0), t > t0, if adj = false

or

   .
  -X(t) = A(t)X(t) + X(t)B(t) + C(t) ,  X(t0) = W(t0), t < t0, if adj = true,

using the periodic generator W determined with the function pgcsylv for the same periodic matrices A, B and C and the same value of the keyword argument adj. The initial time t0 is the nearest time grid value to t, from below, if adj = false, or from above, if adj = true.

The above ODE is solved by employing the integration method specified via the keyword argument solver, together with the required relative accuracy reltol (default: reltol = 1.e-4) and absolute accuracy abstol (default: abstol = 1.e-7). Depending on the desired relative accuracy reltol, lower order solvers are employed for reltol >= 1.e-4, which are generally very efficient, but less accurate. If reltol < 1.e-4, higher order solvers are employed able to cope with high accuracy demands.

The following solvers from the OrdinaryDiffEq.jl package can be selected:

solver = "non-stiff" - use a solver for non-stiff problems (Tsit5() or Vern9());

solver = "stiff" - use a solver for stiff problems (Rodas4() or KenCarp58());

solver = "auto" - use the default solver, which automatically detects stiff problems (AutoTsit5(Rosenbrock23()) or AutoVern9(Rodas5())).

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

Solve the periodic discrete-time Sylvester equation

A'*X*B + isgn*σX = C  for rev = false

or

A*σX*B' + isgn*X = C  for rev = true,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

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)^p, where p is the number of component matrices of X.

Reference:

[1] R. Byers and N. Rhee, “Cyclic Schur and Hessenberg-Schur numerical methods for solving periodic Lyapunov and Sylvester equations,” Dept. of Mathematics, Univ. of Missouri at Kansas City, Tech. Rep., June 1995.

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

pfdsylv(A, B, C; isgn = 1, fast = true) -> X

Solve the forward-time periodic discrete-time Sylvester equation

A*X*B' + isgn*σX = C ,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

prdsylv(A, B, C; isgn = 1, fast = true) -> X

Solve the reverse-time periodic discrete-time Sylvester equation

A'*σX*B + isgn*X = C ,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

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

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

A*X + isgn*σX*B = C  for rev = false

or

isgn*A*σX + X*B = C  for rev = true,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

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.

pfdsylvc(A, B, C; isgn = 1, fast = true) -> X

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

A*X + isgn*σX*B = C ,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.

prdsylvc(A, B, C; isgn = 1, fast = true) -> X

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

isgn*A*σX + X*B = C ,

where σ is the forward shift operator σX(i) = X(i+1) and abs(isgn) = 1.

The periodic matrices A, B and C must have the same type and commensurate periods. The resulting periodic solution X has the period set to the least common commensurate period of A, B and C and the number of subperiods is adjusted accordingly.