Balanced Truncation Model Order Reduction in Limited Frequency and Time Intervals for Discrete-Time Commensurate Fractional-Order Systems

In this paper we investigate an implementation of new model order reduction techniques to linear time-invariant discrete-time commensurate fractional-order state space systems to obtain lower dimensional fractional-order models. Since the models of physical systems correctly approximate the physical phenomena of the modeled systems for restricted time and frequency ranges only, a special attention is given to timeand frequency-limited balanced truncation and frequency-weighted methods. Mathematical formulas for calculation of the timeand frequency-limited, as well as frequency-weighted controllability and observability Gramians, are extended to fractional-order systems. An instructive simulation experiment corroborates the potential of the introduced methodology.


Introduction
In the field of modeling and simulation of fractional-order systems there are two different approaches to application of model order reduction (MOR) techniques: (1) Approximation of fractional-order systems by high integer-order models and their reduction to the low integer-order ones, and (2) reduction of the fractional-order systems without changing the class of the model, i.e., the reduced model is also the fractional-order one.
The first approach can be implemented by either determination of the fractional-order derivative/difference approximators involved in a fractional-order system [1,2] or by selection of integer-order approximators to the whole fractional-order systems [3][4][5][6][7][8]. In both approaches, a very high integer-order model is usually obtained, which is not effective from the computational point of view due to large memory requirements and long simulation times. For these reasons, classical model reduction techniques can be used to reduce integer-order model dimensions [8][9][10]. Therefore, the term MOR for fractional-order systems is usually related to order reduction of integer-order approximators to fractional-order systems.
This paper tackles the issue of MOR for commensurate fractional-order systems where a final result of using the MOR technique is the fractional-order model of lower dimensions. This issue has not been systematically studied and only a few approaches for reduction of continuous-time fractional-order systems have been published [11][12][13][14].
The classical Balanced Truncation (BT) method introduced for classical integer-order systems has been extended to discrete-time fractional-order systems [15]. The reduction paradigms used by the BT method enforce an accurate approximation for the whole range of frequencies. However, models of physical systems, e.g., models for mechanical, electrical and biological systems, characterize a certain adequacy scope, determined by the frequency range, for which the models correctly approximate the physical phenomena of the modeled systems. Likewise, when the reduced model is used to carry out a simulation in the determined time interval, an appropriate approximation accuracy of the output signal y(t) is required only for t lower than a specified final time of simulation. For these reasons, the reduction aims to determine such a reduced model which is particularly accurate in the given frequency range [ω min , ω max ] and/or time interval [t min , t max ]. Such an approach allows for larger errors outside these specified intervals, without negative impact on the usefulness of the obtained reduced model. The time and frequency boundaries can be applied either by using frequency-and time-limited controllability and observability Gramians [16][17][18][19] or by frequency-weighted functions connected to the model which are the subject to the reduction process [20][21][22][23][24][25][26]. In this paper, we focus on the generalization of such approaches to reduction and an accurate approximation in given frequency and time intervals for the discrete-time commensurate fractional-order systems.
The remainder of this paper is structured as follows. A description of fractional-order state space systems considered in the paper is introduced in Section 2. Section 3 includes fundamentals of the MOR concept, in particular for the BT method. Section 4 contains the main result of the paper concerning definitions of controllability and observability Gramians in the time-and frequency-domains for fractional-order systems. Numerical examples of Section 5 illustrate the use of the introduced Gramians in the model reduction process. The paper is completed with the conclusion section.

System Representation
Consider a linear discrete-time commensurate fractional-order (DTCFO) state space system G = {A f , B, C, D} described by where k = 0, 1, . . . , is the discrete time, x(k) ∈ n is the state vector, u(k) and y(k) are the input and output signals, respectively. Matrices A f ∈ n×n , B ∈ n×n u , C ∈ n y ×n , D ∈ n y ×n u describe the system properties, with n u and n y being the numbers of inputs and outputs, respectively. ∆ α x(k + 1) defines the fractional-order difference of order α ∈ (0, 2) which can be represented by the Grünwald-Letnikov fractional difference ( [27], ch. 3.5) with Assuming the zero initial condition, that is ∆ α x(0) = 0, the Z-transform of the system (1) is given by where w(z) is the Z-transform of the "forward-shifted" fractional-order difference where s α is the Laplace transform of the fractional-order derivative and h is the sampling period. Discretized models of continuous-time fractional-order systems where t = kh, can simply be transferred to the system (1) by using the following substitutions Implementation of the Grünwald-Letnikov difference (2) results in increasing computational burden at each time step, which finally becomes computationally infeasible. Therefore, in practical implementations, finite-length expansions are used, for instance finite fractional difference [6,8] with x(l) = 0 ∀ l < 0 and L being the implementation length. It is worth emphasizing that precise approximation of the Grünwald-Letnikov difference with the finite fractional difference needs a very high length L [6].

Model Order Reduction
Let us shift now to the MOR problem for the DTCFO system (1). The fractional-order MOR issue aims towards finding a DTCFO modelG with reduced dimension r < n ∆ αx (k + 1) =Ã fx (k) +Bu(k),x(0) =x 0 y(k) =Cx(k) +Du(k) (9) whereÃ f ∈ r×r ,B ∈ r×n u ,C ∈ n y ×r ,D ∈ n y ×n u ,x(k) ∈ r are such that the approximation errors both in the time domain y(k) −ỹ(k) and in the frequency domain G(z) −G(z) are small for the chosen norm · . In this paper, we concentrate on the BT technique, which is based on the concept of the balanced model realization [29,30], ( [31], ch. 7.1). In order to arrive at the balanced system, the linear state transformation x → Tx is applied to diagonalize the controllability P and observability Q Gramians of the system where σ i , i = 1, ..., n, are the square roots of the eigenvalues for the product of P and Q, which are called the Hankel singular values (HSV). The magnitude of HSV classifies a degree of reachability and observability of states in the system. On this basis, reduction eliminates states corresponding to small values σ i , which means that they have a weak impact on the system properties. The balanced model is obtained by applying transformation matrices in the following way: The order r can be selected on the basis of an approximation error of the reduced system. For the BT method, the H ∞ norm of the approximation error is upper bounded as follows [30], ( [31], ch. 7.2): Calculation of the transformation matrix T is not a unique operation. The exemplary algorithms can be found in References [23,30] ( [31], ch. 7.3), [32].
The BT reduction method enforces an accurate approximation for all times t ∈ (0, ∞) and frequencies ω ∈ . If it is necessary to determine a model (9) which is particularly accurate in a given frequency range [ω min , ω max ] and/or time interval [k min , k max ], then frequency-or time-limited controllability and observability Gramians can be used to calculate the transformation matrix T instead of infinite ones. Higher model accuracy can be obtained now, especially when the optimal values of the parameters of the weighting functions and frequency-/time-intervals are used [33].

Controllability and Observability Gramians for Discrete-Time Fractional-Order Systems
In this section, the definitions of controllability and observability Gramians both in the timeand frequency-domains are recalled. Based on the definitions for integer-order systems, the Gramian generalizations to fractional-order systems are derived here.

Gramians in the Time Domain
The definition of the controllability Gramian is connected to the minimal energy required for the transfer of the system from the zero initial state x(0) = 0 to the final state x(k) = x p , whereas the observability Gramian is defined with relation to the energy generated by the nonzero initial state x(0) = x 0 with the zero input signal u(k) = 0. For asymptotically stable systems, the infinite controllability and observability Gramians are respectively defined as where ξ(k) is the state of the discrete-time system resulting from the input in the form of the Kronecker delta, and η(k) is the output of the system produced by the nonzero initial conditions and the zero input signal. For asymptotically stable integer-order systems, we obtain the well-known formulas for the controllability and observability Gramians, respectively: where A = A f + I. Finally, the Gramians P and Q are the solutions to the discrete-time Lyapunov equations, respectively, are Based on definitions (13), it is easy to formulate the generalized form of the controllability and observability Gramians for the fractional-order systems.

Lemma 1.
Consider an asymptotically stable discrete-time commensurate fractional-order state space system (1), with the Grünwald-Letnikov fractional-order difference (2). Then the controllability and observability Gramians at finite time k L < ∞ are as follows where φ(k), k = 0, 1, ..., are calculated as Proof. The state space equations of the system as in Equation (1) can be rewritten as ( [27], ch. 3.5) The state response for the DTCFO system (17), with the zero initial condition and the Kronecker delta input signal, is as follows: The output response for the DTCFO system (17), with the nonzero initial condition x(0) = x 0 and the zero input signal, is now where The responses (18) and (20) immediately result in Equation (16), which completes the proof.

Remark 2.
As mentioned before, both the minimal energy required for the transfer of the system from the zero initial state to x(k) = x p and the energy generated by the nonzero initial state are obtained for k L → ∞. Implementation of Equation (16) in order to calculate controllability and observability Gramians implies the infinite number of elements φ(k). Furthermore each of φ(k) requires the determination of the Grünwald-Letnikov difference calculated from 0 to k + 1, which is computationally infeasible. Therefore, in contrast to integer-order systems for which the solution of Equation (13) can be determined by solving Lyapunov Equation (15), the Gramians for the fractional-order systems can be calculated for the finite length only.
If the response of the reduced model (9) is expected to match the original fractional-order model (1) in some interval [k 1 , k 2 ], then the balancing of the model can be executed on the basis of the Gramians calculated within this restricted time interval. As the controllability and observability Gramians for the asymptotically stable systems are positive definite, then it can be shown that Therefore, the time-limited Gramians are also positive definite, and within the restricted time interval [k 1 , k 2 ] they can be calculated as follows: where P(k) and Q(k) are given from Equation (16).

Gramians in the Frequency Domain
In addition to the definition (13) given in the time domain, the Gramians can also be expressed in the frequency domain. The connection can be made on the basis of the Plancherel's theorem, which states that the integral of the inner product of two functions in the time domain is equal to the integral of their frequency spectrum. In particular, for the asymptotically stable discrete-time integer-order systems, applying the Plancherel's theorem to Equation (14) Based on definitions (23), it is easy to formulate the generalized form of the controllability and observability Gramians for the fractional-order systems.

Lemma 2.
Consider an asymptotically stable discrete-time commensurate fractional-order state space system (1), with the Grünwald-Letnikov fractional-order difference (2). Then the infinite controllability and observability Gramians of the fractional-order system are respectively given as where w(z) is as in Equation (4), with z = e iθ , θ ∈ [−π, π], and * denotes the complex conjugate transpose.
Proof. For continuous-time fractional-order systems referred to in Remark 1, the input-to-state map becomes (s α I −Ā f ) −1B , while the state-to-output map is C(s α I −Ā f ) −1 . Then the controllability and observability Gramians are as follows: where s = iω, ω ∈ (−∞, ∞). Using the forward-shifted Euler discretization operator as referred to in Remark 1 results immediately in (24), which completes the proof.
If the response of the reduced model is expected to match the full fractional-order model output within a restricted frequency range, then the balancing of the model can be executed on the basis of the Gramians calculated for that specific interval. Similarly, as for time-limited Gramians, the frequency-limited Gramians in restricted frequency interval [Θ 1 , Θ 2 ] can be calculated as follows: where P(Θ) and Q(Θ) are given as

Frequency Weighted Gramians
The required approximation accuracy in the given frequency interval [Θ 1 , Θ 2 ] can also be achieved by implementation of the frequency weighting functions in a form of the external systems connected to the inputs and/or outputs of the full fractional-order model ( [31], ch. 7.6). Such an approach, called Frequency Weighted (FW) method, is a generalization to the BT method designed for asymptotically stable models with asymptotically stable input and output weighting functions with minimal realizations.
Consider an asymptotically stable discrete-time integer-or fractional-order LTI MIMO state space system as an input weighting function and as an output weighting function of orders n i and n o , respectively. Assuming that no pole-zero cancellations occur during the design of the augmented systems GH i and H o G, we arrive at [22,25], ( [31], ch. 7.6) It is well known that the frequency weighted controllability and observability Gramians are computed on the basis of the system connected to the input weight GH i and to the output weight H o G, respectively. (30) consisting of asymptotically stable discrete-time commensurate fractional-order state space system (1), with the Grünwald-Letnikov fractional-order difference (2) and the weighting functions as in Equations (28) and (29). The controllability and observability Gramians of such systems at finite time k L < ∞ are as follows:

Lemma 3. Consider asymptotically stable augmented systems GH i and H o G as in Equation
where φ i (k), φ o (k), k = 0, 1, ..., are calculated in a recurrent way: Proof. The proof directly stems from proof for Lemma 1 with substitution of the system matrices by the matrices of the augmented systems GH i and H o G.
Note that the application of the weighting functions influences the order of the controllability P i and observability Q o Gramians for the augmented systems. Therefore, they are partitioned into two-by-two blocks, so that the dimension of the P 11 ∈ n×n and Q 11 ∈ n×n are the same as the state matrix A f P i = P 11 P 12 P T 12 P 22 , Finally, the frequency-weighted controllabilityP and observabilityQ Gramians of dimensions n × n can be assumed as [20]P = P 11 ,Q = Q 11 (33) The proposed solution, despite its simplicity, may lead to the instability of the reduced model in case of two-sided weighting. Therefore, several modifications to this approach have been proposed to cope with this problem [21][22][23][24][25].
The frequency-weighted Gramians as in Equation (33) can also be defined in the frequency domain. Given that the input-to-state map and the state-to-output map of the fractional-order system are modified by the connected weighting functions, it is possible to generalize the definitions of the infinite controllability and observability Gramians for the fractional-order system (Lemma 2) to the frequency-weighted Gramians. (30), consisting of asymptotically stable discrete-time commensurate fractional-order state space system (1), with the Grünwald-Letnikov fractional-order difference (2) and the weighting functions as in Equations (28) and (29). Then the frequency-weighted controllabilityP and observabilityQ Gramians are defined as

Lemma 4. Consider asymptotically stable augmented systems GH i and H o G as in Equation
Proof. Given that the input-to-state map for the augmented system GH i and the state-to-output map for respectively. TheP andQ are the submatrices of the Gramians for the augmented systems (32) partitioned into two-by-two blocks. Therefore, the proof for the frequency-weighted controllability Gramian follows by noticing that while for the frequency-weighted observability Gramian, Therefore, the frequency-weighted Gramians are the blocks P 11 and Q 11 of the controllability P i and observability Q o Gramians for the augmented systems, respectively, which completes the proof.

Remark 3.
It is important to note that Lemmas 1 to 4 introduce various definitions of controllability and observability Gramians. However, the calculations of the Gramians directly from the above definitions for large-scale DTCFO systems are computationally demanding. In particular, time-domain Gramian definitions as in Lemma 1 and 3 are infeasible for large-scale systems due to the requirement of calculation of the Grünwald-Letnikov difference from 0 to k + 1. The common practice for integer order systems when n u , n y << n is to compute low-rank approximations of the Gramians such that P ≈ SYS T with S ∈ n×l , Y ∈ l×l , l << n. It is motivated by the typical rapid decay of HSV [34,35], which can also be assumed for fractional-order systems. There exists various algorithms for calculating the low-rank Gramian factorizations for integer-order systems [17,19,35,36]. Therefore, future works will be carried out towards the extension of the low-rank approach to increase the efficiency of Gramians calculations for DTCFO systems.

Simulation Examples
In this section, examples of model order reduction for a fractional-order system are presented. All reduced models were obtained by using the BT method with the Gramians calculated within various time-and frequency-intervals, as well as with different frequency-weighted functions. In particular, to calculate transformation matrix T the following Gramians are selected: (1) indefinite Gramians defined in the frequency domain as in Equation (24) -denoted as the BT, (2) frequency-limited Gramians as in Equation (26), denoted as the FLBT, (3) time-limited Gramians as in Equation (16), denoted as the TLBT, (4) frequency-weighted Gramians as in Equation (34), denoted as the FW. The examined discrete-time fractional-order model is an extension of the continuous-time model for a simple mechanical system presented in [16] and (moderate) large-scale dynamical system of order 1006 as in [36].

Example 1.
Consider the DTCFO state space system (6) with the sampling period h = 0.01, fractional order α = 0.85 and Model order reduction is performed from the original six states variables to the reduced four ones. Frequency responses for the full fractional-order system and the reduced models in addition to approximation errors are depicted in Figure 1. Figure 2 shows step responses and their approximation errors for the same models. It is clearly visible that the reduction based on infinite and time-limited Gramians for the time interval t ∈ [0, 10] (s) cannot properly approximate the low-frequency properties of the system. For this reason, in order to improve the approximation for low frequencies, the frequency-interval for frequency-limited Gramians is chosen as Θ ∈ [0, 0.01] (rad/s). For the same purpose, the low-pass Butterworth filter of order n f = 5 and cut-off frequency ω f = 0.01 (rad/s) is selected as a frequency weighted function. Table 1 presents approximation errors for the considered models, where DCE is the steady state approximation error, MSE ω is the mean square approximation error for the frequency responses in the frequency range ω ∈ [10 −3 , 1] (rad/s), H ∞ is the norm approximation error and MSE t is the mean square error for system step response in the discrete-time range t ∈ [0, 100] (s). Figure 3 presents frequency responses and approximation errors for the reduced models obtained in order to improve the quality of approximation for high frequencies. In particular, the third resonance frequency ω = 6.5 (rad/s) is considered. Figure 4 shows the impulse responses and their approximation errors for the same models. For this purpose, the time interval t ∈ [0, 0.1] (s) and frequency interval Θ ∈ [0.1, 100π] (rad/s) are selected for time-and frequency-limited Gramians. Similarly, the frequency weighted function is chosen in a form of a high-pass Butterworth filter of order n f = 5 and cut-off frequency ω f = 0.1 (rad/s). Table 2 presents approximation errors for the analyzed models, in terms of MSE ω for ω > 3 (rad/s), H ∞ -norm and MSE t for system impulse responses within t ∈ [0, 1] (s).        [36] with fractional-order α = 0.95 which is discretized using the sampling period h = 0.002. The calculation of controllability and observability Gramians in the time domain is very computationally demanding for systems of order n = 1006 . Therefore, the reduced models obtained by using only the BT, FLBT and FW methods are compared. All reduced models are of order r = 6. Frequency responses for the full fractional-order system and the reduced models as well as approximation errors are presented in Figure 5. Like in Example 1, it is clearly visible that the reduction based on the infinite Gramians cannot properly approximate the low-frequency properties of the system. For this reason, the frequency-interval for frequency-limited Gramians and frequency weighting functions are chosen the same as in Example 1. In Table 3, approximation errors are listed for the analyzed models, in terms of DCE, MSE ω and H ∞ -norm.  The Matlab scripts used to compute the presented results can be obtained from Supplementary Materials.

Conclusions
This paper presents new results in BT model order reduction in limited time-and frequencyintervals for DTCFO systems. The main contribution of the paper is an introduction of new definitions for controllability and observability Gramians for the fractional-order systems both in the time and frequency domains. These results enable new implementations of the Gramians in the balanced truncation model order reduction method in limited time and frequency intervals as well as in the frequency weighted reduction method. As a result of the reduction process, accurate low-dimension fractional-order approximators in given frequency and/or time intervals can be calculated. Simulation examples confirm the effectiveness of the introduced methodology for order reduction of DTCFO systems.