New integer-order approximations of discrete-time non-commensurate fractional-order systems using the cross Gramian

This paper presents a new approach to modeling of linear time-invariant discrete-time non-commensurate fractional-order single-input single-output state space systems by means of the Balanced Truncation and Frequency Weighted model order reduction methods based on the cross Gramian. These reduction methods are applied to the specific rational (integer-order) FIR-based approximation to the fractional-order system, which enables to introduce simple, analytical formulae for determination of the cross Gramian of the system. This leads to significant decrease of computational burden in the reduction algorithm. As a result, a rational and relatively low-order state space approximator for the fractional-order system is obtained. A simulation experiment illustrates an efficiency of the introduced methodology in terms of high approximation accuracy and low time complexity of the proposed method.

The main problem occurred in fractional-order systems is their infinite time computational complexity, due to the fractional-order derivatives (or differences) incorporated in the model. Therefore, the key step in practical implementations is to approximate fractional-order systems by integer-order models having nearly the same properties in a given frequency range. A number of methods have been proposed in continuous-and discrete-time domains, using transfer function as well as state space representations, e.g., specific finite impulse response (FIR) and infinite impulse response (IIR) representations of fractional-order systems [7,46,[58][59][60]65]. However, an accurate approximation of the fractional-order system in a wide frequency range usually requires a very high integer-order approximator. Thus, to decrease the order of approximators, both in continuous-and discrete-time cases, various model reduction techniques have been proposed [21,36,41,54,60,62].
Among reduction methods, a great attention has been given to the SVD-based methods, which use the balanced model realization theory. This concept enables an easy way of determination of the dominant part of the model and elimination of state variables which have negligible impact on the system properties. In general, the Balanced Truncation (BT) method consists in simultaneous diagonalization of the controllability P and observability Q Gramians [6,23,38,45,60]. A combination of the two, called the cross Gramian, which simultaneously encodes controllability and observability information on the system in a single matrix, can also be used in the model order reduction process. This concept was firstly introduced for single-input single-output (SISO) and symmetric multi-input multi-output (MIMO) systems [17,18,39]. Furthermore, this approach was also extended to nonlinear systems by the use of empirical Gramians [13,24,27].
It is worth mentioning that the reduction process, in particular determination of Gramians for very high order systems, is highly time-consuming. This problem can be solved through the Fourier Model Reduction (FMR) method, which uses discrete-time Fourier coefficients to develop an intermediate order approximation of the original model [49,60,67]. Controllability and observability Gramians can be calculated in an analytical way, which significantly reduces time complexity of the reduction process [49,60]. In this paper, new model order reduction methods combining the FIR-based approximation of a non-commensurate fractional-order system with BT and FW model order reduction methods based on the cross Gramian are proposed.
The remainder of the paper is organized as follows. A representation of noncommensurate fractional-order state space systems and integer-order approximation based on FIR technique are presented in Section 2. Consequently, this Section also gives the background for balancing the systems by use of the cross Gramian and recalls fundamentals of the BT and FW reduction methods for discrete-time state space systems. The main result in terms of new analytical formulae for determination of the cross Gramian for FIR-based approximation of the fractional-order system are presented in Section 3. Numerical examples of Section 4 confirm the effectiveness of the introduced methodology both in terms of modeling accuracy and low time complexity. Conclusions of Section 5 complete the paper.

FIR-based approximation of discrete-time fractional-order system
It is well known that the discrete-time integer-or fractional-order system can be described using the Fourier decomposition with g k ∈ k = 0, 1, . . . , being the impulse response components for the system. In practice, the impulse response components usually approach values very close to zero and can be neglected after a certain implementation length L. In this way, a finite length approximationG L (z) to the G(z) can be obtained. The error bounds corresponding to the introduced approximation usually decide on a value of the chosen implementation length L. However, the error bound given in Ref. [67] is not readily computable. Therefore, one way to determine the appropriate approximation length L is to stop calculating Fourier components when their magnitudes |g k | decrease below a certain value [67]. The way of calculating Fourier components for a commensurate fractional-order system was proposed in Ref. [60] and it can be extended to NCFO systems as follows: where with A, B, C, D and α i as in Eqn. (1) and L being the implementation length. This enables to present the FIR-based approximation of the system for the implementation length L in state space form (G L = {Ã,B,C,D}) as follows: withÃ ∈ L×L ,B ∈ L×1 ,C ∈ 1×L andD ∈ defined as It is worth emphasizing that the state matrixÃ is singular and nilpotent of order L.
Remark 1 For MIMO systems, Fourier components become matrices g k ∈ n y ×n u , where n u and n y denote the numbers of inputs and outputs of the system, respectively. In that case, every element 0 and 1 in the state space form of the FIR-based approximation (6) changes into 0 ∈ n u ×n u and I ∈ n u ×n u being the zero and identity matrices, respectively. Therefore, the order L of the FIR-based MIMO model depends on the number of inputs of the system and implementation length L = Ln u [60,67].

Cross Gramian-based model order reduction
Balancing-related model order reduction methods can be interpreted as performing a linear state transformationx → Tx. This concept is an easy way to determine a dominant part of the model which accurately describes the system dynamics. Transformation (or projection) matrix T and its inverse are not unique. Nevertheless, the majority of algorithms determine its form based on the controllability P and observability Q Gramians of the system [23,29,37,38,49,60,63,66]. Consider a SISO FIR-based approximation of the NCFO system as in Eqn. (6). Let T be the transformation matrix of that system, which diagonalizes the system's controllability and observability Gramians, with decreasing Hankel singular values on the main diagonal where σ 1 ≥ σ 2 ≥ · · · ≥ σ r > σ r+1 ≥ · · · ≥ σ L > 0. Linear state transformation x → Tx enables to calculate the balanced form of the systemḠ L = {Ā,B,C,D} in the following wayḠ where the submatricesÃ r =Ā 11 ∈ r×r ,B r =B 1 ∈ r×1 ,C r =C 1 ∈ 1×r define the reduced-order state space modelG r = {Ã r ,B r ,C r ,D r } of order r < L.
Partitioning T and T −1 in the form of where the matrices T r ∈ r×L , T −1 r ∈ L×r such that T r T −1 r = I r are called truncation matrices and they are the only ones necessary to determine the reduced model [6, ch. 7 The H ∞ -norm of the difference between full (L) and reduced (r) order models is upper bounded by the twice of the sum of the neglected Hankel singular values [6, ch. 7.2] In addition to controllability and observability Gramians, another type of Gramian, known as the cross Gramian X , can be used in the reduction process. The following definitions of reachability R(Ã,B) and observability O(C,Ã) matrices can be a basis for the respective definitions of controllability, observability and cross On this ground, for SISO and symmetric state-space MIMO systems, we can combine the three Gramians' definitions to obtain [18] As it is clearly seen from Eqn. (11), the cross Gramian simultaneously encodes controllability and observability properties into a single matrix X [17,18,39]. For the considered SISO FIR-based approximation (6), due to the nilpotency of the state matrixÃ, the cross Gramian can be expressed as a finite sum and it is always a symmetric matrix for the SISO systems. The controllability and observability Gramians can be calculated from the respective discrete-time Lyapunov equations [ In both cases, square roots of eigenvalues of the product of controllability and observability Gramians and absolute values of eigenvalues of the cross Gramian are input-output invariants and are equal to Hankel singular values of the model For this reason, reduced modelG r can be calculated on the basis of the cross Gramian, making this approach useful in decreasing computational complexity of the reduction algorithm, since solving two Lyapunov equations for P and Q is replaced by solving the single Sylvester equation only.
Various methods can be used to obtain projection matrices from the cross Gramian [4,5,9,27,57]. The first approach can be performed on the basis of singular value decomposition (SVD) of the cross Gramian where U and V are the unitary matrices, U r ∈ L×r , V r ∈ r×L . In the general case, the 1 = diag(σ c 1 , . . . , σ c r ) and 2 = diag(σ c r+1 , . . . , σ c L ) matrices contain σ c i sorted in decreasing order, which are the approximation to the Hankel singular values σ i . In this way, only the approximation to the BT can be obtained [57]. However, for the SISO FIR-based model (6), the cross Gramian is a symmetric matrix, i.e., X = X T (see Eqn. (26)), which provide σ c i = |λ i (X )| = σ i [6, ch. 3.2.3] and for that case the reduction is equivalent to the BT. The truncation matrices T r and T −1 r can be obtained through either the one-sided Galerkin projection or the two-sided Petrov-Galerkin projection The second approach for determination of the transformation matrix T relies on decomposition of matrix X into two blocks, which separate the strongly coupled from the weakly coupled controllable and observable states [3,5,52] where eigenvalues of S l ∈ r×r and S s ∈ (L−r)×(L−r) are the large and small (in magnitude) Hankel singular values of the system. One of the algorithms within the second approach requires of a block-ordered real Schur form of the cross Gramian and calculation of an additional Sylvester equation [3,5,68] where . This enables calculation of the truncation matrices as Another algorithm within the second approach uses orthogonal Givens rotations in order to compute the ordered real Schur forms in ascending and descending orders of absolute eigenvalues along the diagonal [4,52]. Taking into account that the algorithms based on the Schur decomposition are highly computationally involving, they can only be used for moderately complex models. This limits the maximum value of the implementation length L of the FIR-based approximation to the fractional-order model. It finally affects the approximation accuracy of the obtained reduced order model, in particular for low-frequency range. On the other hand, algorithms for computing low-rank approximate solutions to the Sylvester equation, which provide the r-rank approximation to the cross Gramian X , can be used in order to reduce computational complexity and enable approximate balancing reduction of large scale systems [8,9,33,57].

Remark 2
The cross Gramian satisfies Eqn. (11) for SISO and symmetric MIMO systems only. Note that the cross Gramian for non-symmetric, square MIMO systems can be calculated from Eqn. (13), however without any theoretical background as for the symmetric case, so there is no guarantee for obtaining the reduced-order model of appropriate quality [9]. The approaches for non-square MIMO systems rely, e.g., on embedding the system into a symmetric square system by using a symmetrizer matrix [14], which produces a model of the same order but with more inputs and outputs [6, ch. 12.3], [57] or decompose the system into a set of individual SISO subsystems [26,35,43,53].

Frequency weighted model order reduction
The aim of the reduction process based on the above discussed balancing method is system's approximation over all frequencies. In many cases, however, a good approximation is required in designated frequency ranges only. This leads to a specific frequency weighted balanced truncation. The FW method has been developed for stable models and stable input and output weighting functions with minimal realizations W in = {A in , B in , C in , D in } and W out = {A out , B out , C out , D out } of orders m in and m out , respectively. It is well known that, in case of two-sided weighting, the controllability, observability, and cross Gramians are computed based on the system connected with the input weight (G L W in ), the output weight (W outGL ), and both weights (W outGL W in ), respectively [6, ch. 7.6], [16,29,49,63,66,68].
Remark 3 Note that for SISO systems, the input and output weighting functions can be very simply transformed into a single weighting function. Therefore, without the loss of generality for SISO systems, further examination will be conducted for onesided input weighting case only.
Accounting that no pole-zero cancellations occur during forming ofG L W in , the augmented system is given as follows [6, ch. 7.6], [29,49,63,66,68] for which the cross Gramian is the solution of the following Sylvester equatioñ It is important to emphasize that the application of the weighting functions affects the order of the solutions of the Lyapunov/Sylvester equations. However, in order to calculate truncation matrices T r and T −1 r , Gramians of order L are required. A number of algorithms to solve the problem have been developed [16,29,40,63]. For instance, the cross Gramian of the augmented systemG L W in can be partitioned as followsX Like for the Enns technique [16] for frequency weighted controllability and observability Gramians, the frequency weighted cross GramianX can be defined as which means that theX satisfy the following equatioñ A more general approach is based on the Schur complement of a matrix block, which is identical to the method proposed in [40,63] for calculating frequency weighted controllability and observability Gramians where for γ = 0, we have the Enns [16] and for γ = 1 the Lin and Chiu [40] methods. In a general case, when 0 < γ < 1, a combination of these two approaches is chosen [63].

Remark 4
All the subsequent steps of determination of the truncation matrices T r and T −1 r for the FW method are the same as for the methods presented in Section 2.3, with the frequency weighted cross GramianX substituted for X .

Fourier model reduction for continuous time systems
The presented FMR cross Gramian model order reduction method can also be easily applied to reduction of continuous-time fractional-order systems. Note that the FIRbased approximation of the system and reduction algorithm are still realized in the z-domain. Therefore, approximation of continuous-time systems by the use of the considered methodology require the following steps: Taking into account the above items, an additional step in terms of discretization of NCFO system has to be realized.
Consider a stable continuous-time non-commensurate fractional-order LTI SISO with A c ∈ n×n , B c ∈ n×1 , C c ∈ 1×n , D c ∈ denoting the system properties and D α x(t) being the fractional-order derivative vector Fractional-order derivatives D α i x i (t), i = 1, . . . , n can be represented by various definitions involving the Riemann-Liouville, Caputo or Grünwald-Letnikov derivatives [47, ch. 2]. Moreover, to obtain discrete-time non-commensurate fractionalorder system (1), various discretization schemes, such as Euler, Tustin or Al-Alaoui, can be used [2,12]. Taking into account that the simplest discretization procedure can be implemented by the use of the Euler operator to the Grünwald-Letnikov derivative, we arrive at [44, ch. 3.5] where t = kh, k = 0, 1, . . . denotes the (unnormalized) discrete-time and h is the sampling interval. Therefore, transferring from continuous-time system (24) to its discrete-time counterpart (1) results in [44, ch. 3
Taking into account that the reduced discrete-time model (8) is the integer-order model, a conversion process from the discrete to continuous domain can be obtained by using the well known classical methods, e.g., Tustin approach. A value of sampling period h should be chosen with respect to an adequacy range of final reduced order model.

Main result
The main computational cost of the model reduction process based on the cross Gramian is the determination of the Gramian itself. Therefore, the main result of this paper is a new, analytical method for calculation of the cross Gramian of the FIR-based models.

Cross Gramian for FIR-based approximation of NCFO system
In general, the time complexity of calculating the cross Gramian by solving the Sylvester Eq. 13 is a class of cubic one. However, due to the special structure of the FIR-based approximation of discrete-time fractional-order model as in Eqn. (6), it is possible to determine an analytical solution of the cross Gramian.

Lemma 1 Consider the FIR-based approximation (6) of a stable discrete-time noncommensurate fractional-order LTI SISO state space system (1). Then the cross Gramian for that model is an upper triangular Hankel matrix
where g k , k = 1, . . . , L, are as in Eqn. (5).
Proof Because of the special structure of the state space for the FIR-based model, such that the matrixÃ is a lower shift matrix, it results thatÃ k is a k-th diagonal matrix andÃ is nilpotent sinceÃ L = 0. Therefore, the reachability R(Ã,B) and observability O(C,Ã) matrices have finite dimensions (L columns and rows respectively) and are equal to Taking into account (10), we can determine the cross Gramian in a form of Eqn. (26). This completes the proof.

Frequency weighted cross Gramian for FIR-based approximation of NCFO system with input weighting function
The cross Gramian for the augmented systemG L W in is a solution of the Sylvester Eq. 20. However, taking into account (19), the cross Gramian can be determined as a solution to four Sylvester equations of lower dimensions. (30) where X 11 ∈ L×L , X 12 ∈ L×m in , X 21 ∈ m in ×L , X 22 ∈ m in ×m in .

Lemma 2 Consider the discrete-time augmented system as in Eqn. (19). The cross Gramian in a form of Eqn. (21) can be determined as a solution to four consecutive Sylvester equations
Proof Combining Eqns. (19), (20), and (21), we obtain Expanding the matrix description (31) into four submatrix equations, we arrive at Eqns. (27) to (30). This completes the proof.
Due to the special structure ofG L W in given in Eqn. (33), it is possible to present a new, analytical method for determination of the cross Gramian for such a system. (33) consisting of the FIR-based approximation (6) of a stable discrete-time NCFO LTI SISO state space system (1) and the input weighting function (32). Then the cross Gramian for the augmented system is a block matrix such thatX

Theorem 1 Consider the augmented system
where 1) submatrix X 11 ∈ L×L is an upper triangular Hankel matrix with 2) submatrix X 12 ∈ L×m in is as follows 3) submatrix X 21 ∈ m in ×L is a rectangular Hankel matrix . . .
4) submatrix X 22 ∈ m in ×m in is as follows and, for all the four items above, g k = 0 for k > L, w k = 0 for k > m in ,r k = 0 for k < 1 or k > L.
Proof Consider the cross Gramian of the augmented systemG L W in partitioned as in Eqn. (34), where the submatrices ofX can be concisely presented as follows: with the dimensions β = γ = L for X 11 , β = L, γ = m in for X 12 , β = m in , γ = L for X 21 and β = γ = m in for X 22 resulting from Lemma 2.
SinceÃ and A in are lower shift matrices (compare Eqn. (33)), the left hand sides of the four consecutive Sylvester Eqs. 27 to (30) are as follows: The right hand sides of the Sylvester Eqs. 27 and (28) are respectively equal to Both above matrices contain nonzero entries in the first rows only, which proves that X 11 and X 21 are the upper triangular Hankel matrix and rectangular Hankel matrix, respectively.
Accounting for the elements x 12 i,j and x 22 i,j occurring in the previous rows, e.g., starting from x 12 1,2 = σ 2 in X 12 and from x 22 1,2 = g 0 w 2 in X 22 , we arrive at Eqns. (36) and (38). This completes the proof.
As mentioned above, the cross Gramian for the augmented system, consisting of the FIR-based approximation for NCFO system and input weighting function, satisfies (20) for SISO systems only. Note that the cross Gramian for non-symmetric square MIMO systems can be determined using Theorem 1 with substitution of all the elements g k , w k by the respective matrices g k ∈ n u ×n u , w k ∈ n u ×n u and all the elements 0 and 1 by the respective zero and identity matrices 0 ∈ n u ×n u and I ∈ n u ×n u , where n u denotes both the number of inputs and outputs of the system. However, for non-symmetric square MIMO systems the singular values of the cross Gramian σ c i (14) are only the approximation to the Hankel singular values σ i [57]. Therefore, it must be emphasized that without any theoretical background as for the SISO case, there is no guarantee for obtaining the reduced-order model of appropriate quality.

Frequency weighted Fourier model reduction algorithm
The frequency weighted Fourier model reduction (FMR-FW) algorithm for reduction of fractional-order systems can be summarized in the following steps: 1. For continuous-time NCFO system, use discretization procedure to obtain discrete-time counterpart as in Eqn. (25). 2. Calculate g k components for NCFO system using Eqn. (5). 3. Specify w k components of the input weighting function (32). 4. Create the augmented systemG L W in as in Eqn. (33), consisting of the FIR-based approximation of NFCO system and input weighting function. 5. Use Theorem 1 to analytically calculate the cross GramianX of the system G L W in . 6. Select coefficient γ and calculate frequency weighted cross GramianX based on Eqn. (23). 7. Calculate the truncation matrices T r and T −1 r using either the one-sided Galerkin projection (15) or the two-sided Petrov-Galerkin projection (16) or the algorithm based on the Schur decomposition (18). 8. Reduce order of the model as in Eqn. (8). 9. UseG r as a final model for the discrete-time NCFO system or convertG r to the continuous-time domain using classical methods (e.g., Tustin approach) in order to obtain a final model for the continuous-time NCFO system.

Simulation examples
The effectiveness of the reduction process by the use of the introduced methodology will be demonstrated on test examples. All the computations were done on a board 2 IntelXeon E5-2683V3 CPUs with 2.0 GHz clock and 128 GB RAM using MATLAB/Simulink 2015b software.
Example 1 Consider the discrete-time NCFO state space system (1) with The system is approximated by the FIR-based models (6), with the H ∞ -norm approximation errors equal to 0.0855, 0.0204, 0.0044 for the respective implementation lengths L being 100, 1000, and 10000. Hankel singular values for all models are presented in Fig. 1a. Figure 1b shows approximation errors, the upper bound of H ∞norm (9) calculated on the basis of HSV for FIR-based models and the actual value of H ∞ -norm for the reduced models obtained by use of the FMR-BT method. The results show that approximation accuracy for the FIR-based model plays a crucial role in the reduction process as the global approximation error for the reduced model is a sum of errors for the FIR-based approximation and the reduction process.
Example 2 Consider the NFCO system as in Example 1 approximated by the FIRbased model (6), with the implementation length L = 10 4 . The input frequency weighting function is applied in the form of a low-pass FIR digital filter with cutoff frequency ω c = 0.01π [rad/s] and m in = 12. The system is reduced to models of order r = 4 using the FMR-BT method based on the cross Gramian given as in Lemma 1 as well as the FMR-FW method based on the frequency weighted cross Gramian as in Eqn. (23) for γ = 0, 0.35, 0.7, 1.0 and the one-sided Galerkin projection for determination of the truncation matrices T r and T −1 r . The cross Gramian of the augmented systemX is calculated based on Theorem 1. Since the matrix X 22 may be ill-conditioned, the Moore-Penrose pseudoinverse can be used. The frequency responses of NCFO system, high integer-order FIR-based approximation and reduced models, as well as the approximation errors are presented in Fig. 2.
Example 3 Consider the same example of fractional-order system and the FMR-FW algorithm as in Example 2. Figure 3 shows the approximation errors of reduced models obtained by the use two algorithms of determination for truncation matrices T r and T −1 r , that is (1) the one-sided Galerkin projection (15), denoted as FMR-FW-G and (2) the algorithm based on the Schur decomposition (18), denoted as FMR-FW-S. Table 1 presents the values of the approximation errors for the analyzed models in terms of DCE -steady-state approximation error, MSE -mean square approximation error for the frequency characteristics in the frequency range ω ∈ [10 −5 , π] and H ∞ -norm approximation error.  To compare time complexities of (each iteration of) the FMR-FW algorithm the fractional-order system was reduced for various implementation lengths L and various orders m in of the input weighting function. Table 2 shows times required for determination of the augmented system (33), item 1 in the table, calculation of the cross Gramian from the Sylvester Eq. 20, item 2, calculation of the cross Gramian based on Lemma 2, item 3, calculation of the cross Gramian based on Theorem 1, item 4, calculation of truncation matrices T r and T −1 r based on the one-sided Galerkin projection (15), item 5, calculation of truncation matrices based on the real Schur decomposition (18), item 6.
The results presented in Figs. 2 and 3 and Table 1 show that FMR-FW gives a good approximation accuracy, in particular for low frequency range, in comparison to the FMR-BT method, due to application of the low-pass filter. It can be seen that the approximation error in the low frequency range is in a large part the result of the FIR-based approximation of the NCFO system. Reduced order models obtained by the use of both methods for determination of truncation matrices give very similar results. Taking into account that the algorithm based on the Schur decomposition is much more computationally involving, it is strongly recommended to use the onesided Galerkin projection. It can also be seen from Table 2 that the computational times for the FMR-FW algorithm using the introduced analytical solution to the cross Gramian and the one-sided Galerkin projection is the least time consuming.
The Matlab scripts used to compute the presented results can be obtained from: http://doi.org/10.5281/zenodo.1256699.

Conclusions
In this paper, a new approximation methodology for modeling of discrete-time non-commensurate fractional-order state space systems by integer-order state-space models has been offered. The introduced FW reduction method is based on cross Gramian of a high-order FIR-based approximation to the fractional-order system. The main contribution of the paper is the introduction of analytical formulae enabling low-cost calculation of the cross Gramians, including the frequency weighted one, instead of a very time consuming process of classical solving the Sylvester equation. Simulation examples illustrate the efficiency of the proposed methodology, both in terms of high accuracy within prespecified frequency ranges and low computation effort of the reduction algorithms.