Periodic Riccati equation solvers

pcric(A, R, Q; K = 10, N = 1, adj = false, solver, reltol, abstol, fast, intpol, dt) -> (X, EVALS)

Solve the periodic Riccati differential equation

.                                                 
X(t) = A(t)X(t) + X(t)A(t)' + Q(t) - X(t)R(t)X(t) ,  if adj = false,

or

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

and compute the stable closed-loop characteristic multipliers in EVALS (see pgcric for details).

The periodic matrices A, R and Q must have the same type, the same dimensions and commensurate periods, and additionally R and Q 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, R and Q and the number of subperiods is adjusted accordingly. Note: Presently the PeriodicTimeSeriesMatrix and PeriodicSwitchingMatrix types are not supported.

If fast = true (default) the multiple-shooting method is used in conjunction with fast pencil reduction techniques, as proposed in [1], to determine the periodic solution in t = 0 and a multiple point generator of the appropriate periodic differential Riccati equation is determined (see [2] for details). If fast = false, the multiple-shooting method is used in conjunction with the periodic Schur decomposition to determine multiple point generators directly from the stable periodic invariant subspace of an appropriate symplectic transition matrix (see also [2] for more details).

For the determination of the multiple point periodic generators an implicit Runge-Kutta Gauss-Legendre 16th order method from the IRKGaussLegendre.jl package is employed to integrate the appropriate Hamiltonian system [2].

For the evaluation of solution via interpolation or ODE integration, the following solvers from the OrdinaryDiffEq.jl package can be selected using the keyword argument solver:

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

The accuracy of the computed solutions can be controlled via the relative accuracy keyword reltol (default: reltol = 1.e-4) and absolute accuracy keyword 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.

For large values of K, parallel computation can be used to determine the matrices of the discrete-time problem or the multiple domain interpolation based solutions. This requires to start 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. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

pgcric(A, R, Q[, K = 1]; adj = false, solver, reltol, abstol, fast, PSD_SLICOT) -> (X, EVALS)

Compute periodic generators for the periodic Riccati differential equation in the filtering form

.                                                  
X(t) = A(t)X(t) + X(t)A(t)' + Q(t) - X(t)R(t)X(t), if adj = false,

or in the control form

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

where A(t), R(t) and Q(t) are periodic matrices of commensurate periods, with A(t) square, R(t) symmetric and positive definite, and Q(t) symmetric and positive semidefinite. 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(t), R(t) and Q(t) are constant matrices and N = K otherwise. EVALS contains the stable characteristic multipliers of the monodromy matrix of the corresponding Hamiltonian matrix (also called closed-loop characteristic multipliers). The period of X is set to the least common commensurate period of A(t), R(t) and Q(t) 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.

If fast = true (default) the multiple-shooting method is used in conjunction with fast pencil reduction techniques, as proposed in [1], to determine the periodic solution in t = 0 and a multiple point generator of the appropriate periodic differential Riccati equation is determined (see [2] for details). If fast = false, the multiple-shooting method is used in conjunction with the periodic Schur decomposition to determine multiple point generators directly from the stable periodic invariant subspace of an appropriate symplectic transition matrix (see also [2] for more details).

The recommended solver to integrate the appropriate Hamiltonian system [2] is the implicit Runge-Kutta Gauss-Legendre 16th order method from the IRKGaussLegendre.jl package, which can be selected with solver = "symplectic" (default).

Other ODE solvers from the OrdinaryDiffEq.jl package can be also 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 = "linear" - use a special solver for linear ODEs (MagnusGL6()) with fixed time step dt;

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

The accuracy of the computed solutions can be controlled via the relative accuracy keyword reltol (default: reltol = 1.e-4) and absolute accuracy keyword 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.

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. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

tvcric_eval(t, W, A, R, Q; adj, solver, reltol, abstol, dt) -> Xval

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

  .                                                
  X(t) = A(t)X(t) + X(t)A(t)' + Q(t) - X(t)R(t)X(t) ,  X(t0) = W(t0), t0 < t, if adj = false (default),

or

  .                                                
 -X(t) = A(t)'X(t) + X(t)A(t) + Q(t) - X(t)R(t)X(t) ,  X(t0) = W(t0), t0 > t, if adj = true,

using the periodic generator W determined with the function pgcric for the same periodic matrices A, R and Q 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 resulting Xval is a symmetric matrix.

The ODE solver to be employed can be specified using the keyword argument solver, (default: solver = "non-stiff") 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())).

prcric(A, B, R, Q; K = 10, N = 1, solver, intpol, reltol, abstol, fast) -> (X, EVALS, F)

Compute the symmetric stabilizing solution X(t) of the periodic control related Riccati differential equation

 .                                                -1 
-X(t) = A(t)'X(t) + X(t)A(t) + Q(t) - X(t)B(t)R(t)  B(t)'X(t) ,

the periodic stabilizing state-feedback gain

           -1
F(t) = R(t)  B(t)'X(t)

and the corresponding stable characteristic multipliers EVALS of A(t)-B(t)F(t).

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

If fast = true (default) the multiple-shooting method is used in conjunction with fast pencil reduction techniques, as proposed in [1], to determine the periodic solution in t = 0 and a multiple point generator of the appropriate periodic differential Riccati equation is determined (see [2] for details). If fast = false, the multiple-shooting method is used in conjunction with the periodic Schur decomposition to determine multiple point generators directly from the stable periodic invariant subspace of an appropriate symplectic transition matrix (see also [2] for more details).

The keyword argument K specifies the number of grid points to be used for the resulting multiple point periodic generator (default: K = 10). 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 interplation-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.

For the determination of the multiple point periodic generators an implicit Runge-Kutta Gauss-Legendre 16th order method from the IRKGaussLegendre.jl package is employed to integrate the appropriate Hamiltonian system [2].

For the evaluation of solution via interpolation or ODE integration, the following solvers from the OrdinaryDiffEq.jl package can be selected using the keyword argument solver:

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

The accuracy of the computed solutions can be controlled via the relative accuracy keyword reltol (default: reltol = 1.e-4) and absolute accuracy keyword 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.

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. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

pfcric(A, C, R, Q; K = 10, N = 1, solver, intpol, reltol, abstol, fast) -> (X, EVALS, F)

Compute the symmetric stabilizing solution X(t) of the periodic filtering related Riccati differential equation

.                                                 -1 
X(t) = A(t)X(t) + X(t)A(t)' + Q(t) - X(t)C(t)'R(t)  C(t)X(t) ,

the periodic stabilizing Kalman gain

                    -1
F(t) = X(t)C(t)'R(t)

and the corresponding stable characteristic multipliers EVALS of A(t)-F(t)C(t).

The periodic matrices A, C, R and Q must have the same type and commensurate periods, and additionally R must be symmetric positive definite and Q must be symmetric positive semidefinite. 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, C, R and Q and the number of subperiods is adjusted accordingly.

If fast = true (default) the multiple-shooting method is used in conjunction with fast pencil reduction techniques, as proposed in [1], to determine the periodic solution in t = 0 and a multiple point generator of the appropriate periodic differential Riccati equation is determined (see [2] for details). If fast = false, the multiple-shooting method is used in conjunction with the periodic Schur decomposition to determine multiple point generators directly from the stable periodic invariant subspace of an appropriate symplectic transition matrix (see also [2] for more details).

The keyword argument K specifies the number of grid points to be used for the resulting multiple point periodic generator (default: K = 10). 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 interplation-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.

For the determination of the multiple point periodic generators an implicit Runge-Kutta Gauss-Legendre 16th order method from the IRKGaussLegendre.jl package is employed to integrate the appropriate Hamiltonian system [2].

For the evaluation of solution via interpolation or ODE integration, the following solvers from the OrdinaryDiffEq.jl package can be selected using the keyword argument solver:

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

The accuracy of the computed solutions can be controlled via the relative accuracy keyword reltol (default: reltol = 1.e-4) and absolute accuracy keyword 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.

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. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

 prdric(A, B, R, Q[, S]; itmax = 0, nodeflate = false, fast, rtol) -> (X, EVALS, F)

Solve the periodic Riccati difference equation

  X(i) = Q(i) + A(i)'X(i+1)A(i) - (A(i)'X(i+1)B(i) + S(i))*
                                 -1
         (B(i)'X(i+1)B(i) + R(i))  (A(i)'X(i+1)B(i) + S(i))'

and compute the stabilizing periodic state feedback

                                  -1
  F(i) = -(B(i)'X(i+1)B(i) + R(i))  (B(i)'X(i+1)A(i) + S(i)')

and the corresponding stable closed-loop characteristic multipliers of A(i)-B(i)F(i) in EVALS.

The n×n and n×m periodic matrices A(i) and B(i) are contained in the PeriodicArray objects A and B, and must have the same sampling time. R(i), Q(i) and S(i) are m×m, n×n and n×m periodic matrices of same sampling times as A and B, and such that R(i) and Q(i) are symmetric. R(i), Q(i) and S(i) are contained in the PeriodicArray objects R, Q and S. R, Q and S can be alternatively provided as constant real matrices. The resulting symmetric periodic solution X and periodic state feedback gain F have the period set to the least common commensurate period of A, B, R and Q and the number of subperiods is adjusted accordingly.

If fast = true, the fast structure exploiting pencil reduction based method of [1] is used to determine a periodic generator in X(1), which allows to generate iteratively the solution over the whole period. If fast = false (default), the periodic Schur decomposition based approach of [1] is employed, applied to a symplectic pair of periodic matrices. If nodeflate = false (default), the underlying periodic pencil is preprocessed to eliminate (deflate) the infinite characteristic multipliers originating from the problem structure. If nodeflate = true, no preliminary deflation is performed.

An iterative refining of the accuracy of the computed solution can be performed by using itmax = k, with k > 0 (default: k = 0).

To detect singularities of involved matrices, the keyword parameter rtol = tol can be used to specify the lower bound for the 1-norm reciprocal condition number. The default value of tol is n*ϵ, where ϵ is the working machine epsilon.

References

[1] A. Varga. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

 prdric(A, B, R, Q[, S]; itmax = 0, nodeflate = false, fast, rtol) -> (X, EVALS, F)

Solve the periodic Riccati difference equation

  X(i) = Q(i) + A(i)'X(i+1)A(i) - (A(i)'X(i+1)B(i) + S(i))*
                                 -1
         (B(i)'X(i+1)B(i) + R(i))  (A(i)'X(i+1)B(i) + S(i))'

and compute the stabilizing periodic state feedback

                                  -1
  F(i) = -(B(i)'X(i+1)B(i) + R(i))  (B(i)'X(i+1)A(i) + S(i)')

and the corresponding stable closed-loop core characteristic multipliers of A(i)-B(i)F(i) in EVALS.

The n(i+1)×n(i) and n(i+1)×m periodic matrices A(i) and B(i) are contained in the PeriodicMatrix objects A and B, and must have the same sampling time. R(i), Q(i) and S(i) are m×m, n(i)×n(i) and n(i)×m periodic matrices of same sampling times as A and B, and such that R(i) and Q(i) are symmetric. R(i), Q(i) and S(i) are contained in the PeriodicMatrix objects R, Q and S. R, Q and S can be alternatively provided as constant real matrices. The resulting n(i)×n(i) symmetric periodic solution X(i) and m×n(i) periodic state feedback gain F(i) have the period set to the least common commensurate period of A, B, R and Q and the number of subperiods is adjusted accordingly.

If fast = true, the fast structure exploiting pencil reduction based method of [1] is used to determine a periodic generator in X(j), which allows to generate iteratively the solution over the whole period. The value of j corresponds to the least dimension nc of n(i) (which is also the number of core characteristic multipliers). If fast = false (default), the periodic Schur decomposition based approach of [1] is employed, applied to a symplectic pair of periodic matrices. If nodeflate = false (default), the underlying periodic pencil is preprocessed to eliminate (deflate) the infinite characteristic multipliers originating from the problem structure. If nodeflate = true, no preliminary deflation is performed.

An iterative refining of the accuracy of the computed solution can be performed by using itmax = k, with k > 0 (default: k = 0).

To detect singularities of involved matrices, the keyword parameter rtol = tol can be used to specify the lower bound for the 1-norm reciprocal condition number. The default value of tol is n*ϵ, where ϵ is the working machine epsilon and n is the maximum of n(i).

References

[1] A. Varga. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

 pfdric(A, C, R, Q[, S]; itmax = 0, nodeflate = false, fast, rtol) -> (X, EVALS, F)

Solve the periodic Riccati difference equation

  X(i+1) = Q(i) + A(i)X(i)A(i)' - (A(i)X(i)C(i)' + S(i))*
                                 -1
           (C(i)X(i)C(i)' + R(i))  (A(i)X(i)C(i)' + S(i))'

and compute the stabilizing periodic Kalman gain

                                                       -1
  F(i) = -(C(i)X(i)A(i)' + S(i)')(C(i)X(i)C(i)' + R(i))

and the corresponding stable Kalman filter characteristic multipliers of A(i)-F(i)C(i) in EVALS.

The n×n and m×n periodic matrices A(i) and C(i) are contained in the PeriodicArray objects A and C, and must have the same sampling time. R(i), Q(i) and S(i) are m×m, n×n and m×n periodic matrices of same sampling times as A and C, and such that R(i) and Q(i) are symmetric. R, Q and S can be alternatively provided as constant real matrices. The resulting symmetric periodic solution X and Kalman filter gain F have the period set to the least common commensurate period of A, C, R and Q and the number of subperiods is adjusted accordingly.

The dual method of [1] is employed (see prdric for the description of keyword parameters).

References

[1] A. Varga. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.

 pfdric(A, C, R, Q[, S]; itmax = 0, nodeflate = false, fast, rtol) -> (X, EVALS, F)

Solve the periodic Riccati difference equation

  X(i+1) = Q(i) + A(i)X(i)A(i)' - (A(i)X(i)C(i)' + S(i))*
                                 -1
           (C(i)X(i)C(i)' + R(i))  (A(i)X(i)C(i)' + S(i))'

and compute the stabilizing periodic Kalman gain

                                                       -1
  F(i) = -(C(i)X(i)A(i)' + S(i)')(C(i)X(i)C(i)' + R(i))

and the corresponding stable Kalman filter core characteristic multipliers of A(i)-F(i)C(i) in EVALS.

The n(i+1)×n(i) and m×n(i) periodic matrices A(i) and C(i) are contained in the PeriodicMatrix objects A and C, and must have the same sampling time. R(i), Q(i) and S(i) are m×m, n(i)×n(i) and m×n(i) periodic matrices of same sampling times as A and C, and such that R(i) and Q(i) are symmetric. R, Q and S can be alternatively provided as constant real matrices. The resulting symmetric periodic solution X and Kalman filter gain F have the period set to the least common commensurate period of A, C, R and Q and the number of subperiods is adjusted accordingly.

The dual method of [1] is employed (see prdric for the description of keyword parameters).

References

[1] A. Varga. On solving periodic Riccati equations. Numerical Linear Algebra with Applications, 15:809-835, 2008.