Periodic Lyapunov equation solvers

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

Solve the periodic Lyapunov differential equation

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

or

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

The periodic matrices A and C must have the same type, the same dimensions and commensurate periods. Additionally C must be symmetric. The resulting symmetric 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 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 Lyapunov equation into a discrete-time periodic Lyapunov 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 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 Lyapunov 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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

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

Solve the periodic reverse-time Lyapunov differential equation

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

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

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

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

Solve the periodic forward-time Lyapunov differential equation

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

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

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

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

Compute periodic generators for the periodic Lyapunov differential equation

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

or

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

The periodic matrices A and C must have the same type, the same dimensions and commensurate periods, and additionally C must be symmetric. K specifies the number of grid points to be used for the discretization of the continuous-time problem (default: K = 1). If A 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 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 and C are constant matrices and N = K otherwise. If A 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 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 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 Lyapunov differential equation into a discrete-time periodic Lyapunov equation satisfied by the generator solution in the grid points and then to compute the solution by solving an appropriate discrete-time periodic Lyapunov 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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

pgclyap2(A, C, E, [, K = 1]; solver, reltol, abstol) -> (X,Y)

Compute the solutions of the periodic differential Lyapunov equations

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

and

 .
-Y(t) = A(t)'*Y(t) + Y(t)*A(t) + E(t).

The periodic matrices A, C and E must have the same dimensions, the same type and commensurate periods. Additionally C and E must be symmetric. The allowed types are PeriodicFunctionMatrix, PeriodicSymbolicMatrix and HarmonicArray.

K specifies the number of grid points to be used for the discretization of the continuous-time problem (default: K = 1). If A, C and E have the types PeriodicFunctionMatrix, HarmonicArray, FourierFunctionMatrix or PeriodicSymbolicMatrix, then the resulting X and Y are collections of periodic generator matrices determined as periodic time-series matrices with N components, where N = 1 if A, C and E are constant matrices and N = K otherwise. The period of X and Y is set to the least common commensurate period of A, C and E and the number of subperiods is adjusted accordingly. Any component matrix of X or Y 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 Lyapunov equations into discrete-time periodic Lyapunov equations satisfied by the generator solutions in the grid points and then to compute the solutions by solving appropriate discrete-time periodic Lyapunov equations using the periodic Schur method of [2].

The ODE solver to be employed to convert the continuous-time problems into discrete-time problems 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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

pgclyap2(A, C::AbstractMatrix, E, [, K = 1]; solver, reltol, abstol) -> (X,Y)

Compute the solution of the discrete-time periodic Lyapunov equation

X(i+1) = Φ(i)*X(i)*Φ'(i) + W(i), i = 1, ..., K, X(K+1) := X(1)

and a periodic generator for the periodic Lyapunov differential equations

 .
-Y(t) = A(t)'Y(t) + Y(t)A(t) + E(t).

The periodic matrices A and E and the constant matrix C must have the same dimensions, and A and E must have the same type and commensurate periods. Additionally C and E must be symmetric. Φ(i) denotes the transition matrix on the time interval [Δ*(i-1), Δ*i] corresponding to A, where Δ = T/K with T the common period of A and E. W(i) = 0 for i = 1, ..., K-1 and W(K) = C. If A and E have the types PeriodicFunctionMatrix, HarmonicArray, FourierFunctionMatrix or PeriodicSymbolicMatrix, then the resulting Y is a collection of periodic generator matrices determined as a periodic time-series matrix with N components, where N = 1 if A and E are constant matrices and N = K otherwise. The period T of Y is set to the least common commensurate period of A and E and the number of subperiods is adjusted accordingly. Any component matrix of Y 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 Lyapunov into a discrete-time periodic Lyapunov equation satisfied by the generator solution in the grid points and then to compute the solution by solving an appropriate discrete-time periodic Lyapunov equation using the periodic Schur method of [2].

The ODE solver to be employed to convert the continuous-time problems into discrete-time problems 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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

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

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

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

or

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

using the periodic generator W determined with the function pgclyap for the same periodic matrices A 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())).

pcplyap(A, C; K = 10, adj = false, solver, reltol, abstol) -> U

Compute the upper triangular periodic factor U(t) of the solution X(t) = U(t)U(t)' of the periodic Lyapunov differential equation

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

or of the solution X(t) = U(t)'U(t) of the periodic Lyapunov differential equation

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

The periodic matrices A and C must have the same type, commensurate periods and A must be stable. The resulting upper triangular periodic factor U has the type PeriodicFunctionMatrix and U(t) can be used to evaluate the value of U at time t. U has the period set to the least common commensurate period of A and C and the number of subperiods is adjusted accordingly.

An extension of the multiple-shooting method of [1] is employed to convert the (continuous-time) periodic differential Lyapunov equation into a discrete-time periodic Lyapunov equation satisfied by a multiple point 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). The upper triangular factor of the multiple point generator is computed by solving the appropriate discrete-time periodic Lyapunov equation using the iterative method (Algorithm 5) of [2]. The resulting periodic generator is finally converted into a periodic function matrix which determines for a given t the function value U(t) by integrating the appropriate ODE from the nearest grid point value.

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

To speedup function evaluations, interpolation based function evaluations can be used by setting the keyword argument intpol = true (default: intpol = true). Interpolation is not possible if A and C are of type PeriodicSwitchingMatrix.

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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

prcplyap(A, C; K = 10, solver, reltol, abstol) -> U

Compute the upper triangular periodic factor U(t) of the solution X(t) = U(t)'U(t)

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

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

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

pfcplyap(A, B; K = 10, solver, reltol, abstol) -> U

Compute the upper triangular periodic factor U(t) of the solution X(t) = U(t)U(t)'

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

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

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

pgcplyap(A, C[, K = 1]; adj = false, solver, reltol, abstol) -> U

Compute upper triangular periodic generators U(t) of the solution X(t) = U(t)U(t)' of the periodic Lyapunov differential equation

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

or of the solution X(t) = U(t)'U(t) of the periodic Lyapunov differential equation

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

The periodic matrices A and C must have the same type, commensurate periods and A must be stable. The resulting U is a collection of periodic generator matrices determined as a periodic time-series matrix with N components, where N = 1 if A and C are constant matrices and N = K otherwise. The period of U is set to the least common commensurate period of A and C and the number of subperiods is adjusted accordingly. Any component matrix of U is a valid initial value to be used to generate the solution over a full period by integrating the appropriate differential equation. An extension of the multiple-shooting method of [1] is employed, first, to convert the continuous-time periodic Lyapunov into a discrete-time periodic Lyapunov equation satisfied by the generator solution in K time grid points and then to compute the solution by solving an appropriate discrete-time periodic Lyapunov equation using the iterative method (Algorithm 5) 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())).

For large values of K, 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] A. Varga. Periodic Lyapunov equations: some applications and new algorithms. Int. J. Control, vol, 67, pp, 69-87, 1997.

 tvcplyap_eval(t, U, A, C; adj = false, solver, reltol, abstol) -> Uval

Compute the time value Uval := U(t) of the upper triangular periodic generators U(t) of the solution X(t) = U(t)U(t)' of the periodic Lyapunov differential equation

  .
  X(t) = A(t)X(t) + X(t)A(t)' + C(t)C(t)' , X(t0) = U(t0)U(t0)', t > t0, if adj = false,

or of the solution X(t) = U(t)'U(t) of the periodic Lyapunov differential equation

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

using the periodic generator U determined with the function pgcplyap for the same periodic matrices A 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), absolute accuracy abstol (default: abstol = 1.e-7) and stepsize dt (default: dt = 0, only used if solver = "symplectic"). 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())).

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

Solve the periodic discrete-time Lyapunov equation

A'σXA + C = X  for adj = true

or

AXA' + C = σX  for adj = false,

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

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

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.

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

Solve the pair of periodic discrete-time Lyapunov equations

AXA' + C  = σX, 
A'σYA + E = Y,

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

The periodic matrices A, C and E must have the same type, the same dimensions and commensurate periods, and additionally C and E must be symmetric. The resulting symmetric periodic solutions X and Y have the period set to the least common commensurate period of A, C and E and the number of subperiods is adjusted accordingly.

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.

prdlyap(A, C; stability_check = false) -> X

Solve the reverse-time periodic discrete-time Lyapunov equation

A'σXA + C = X

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

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

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.

pfdlyap(A, C; stability_check = false) -> X

Solve the forward-time periodic discrete-time Lyapunov equation

AXA' + C = σX

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

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

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.

pdplyap(A, C; adj = true) -> 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 must have the same type, the same dimensions and commensurate periods. The resulting upper triangular periodic matrix U has the period set to the least common commensurate period of A and C and the number of subperiods is adjusted accordingly.

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.

prdplyap(A, C) -> U

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

A'σXA + C'C = X

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 must have the same type, the same dimensions and commensurate periods. The resulting upper triangular periodic matrix U has the period set to the least common commensurate period of A and C and the number of subperiods is adjusted accordingly.

Note: X is the observability Gramian of the periodic pair (A,C).

pfdplyap(A, B) -> U

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

AXA' + BB' = σX

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 B must have the same type, the same dimensions and commensurate periods. The resulting upper triangular periodic matrix U has the period set to the least common commensurate period of A and B and the number of subperiods is adjusted accordingly.

Note: X is the reachability Gramian of the periodic pair (A,B).