A Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems

This paper presents a robust model order reduction technique with guaranteed stability, minimum phase, and matched steady-state response for linear time-invariant single-input-single-output systems. The proposed approach is generalized, allowing the designer to select any desired order of the reduced-order model (ROM). In contrast to the published literature, which primarily uses the time-domain behavior, the proposed technique utilizes the frequency-domain information of the full-order system. The suggested strategy allows the determination of the optimal ROM in a single step, simpler than the various recently reported mixed methods. The robustness is demonstrated using convergence studies and statistical measures about the final solution quality and model coefficients. The superiority over the recent literature is illustrated through four numerical examples using various time-domain and frequency response performance metrics.


I. INTRODUCTION
M ODEL order reduction (MOR) deals with the approximation of complex, high-order models with the loworder ones [1]. A reduced-order model (ROM) requires fewer variables and parameters for representation either in transfer function or state-space form while preserving the essential characteristics of the original system.
The transfer function of a ROM can be defined by (2): where A k (k = 0, 1, . . . , M) and B k (k = 0, 1, . . . , N−1) are the coefficients of the numerator and denominator polynomials, respectively, of G R (s); M and N are integers; M < m, and N < n.
ROM helps alleviate the computational burden during process control simulation studies, reduces simulation time during numerical analysis of dynamical systems, and reduces hardware overhead for practical implementation. Research in the field of MOR has remained vigorously active in the past decade with applications reported in various fields of science and engineering, such as semiconductor components packaging [2], electromagnetism [3], etc.
MOR strategies for LTI SISO systems, such as the step and impulse inputs based response matching [4], mixed positive-bounded balanced truncation [5], and the balanced realization technique (BRT) [6], have been reported in the literature. Preservation of ROM passivity was guaranteed using the truncated BRT algorithms [7]. The pole clustering technique was employed in [8], while clustering along with the retention of dominant poles was reported in [9]. Geo-metric programming [10], the least-squares method [11], and moment matching techniques [12], [13] have been explored to solve MOR problems. Hankel matrices were used for the reduction of large-scale models in the time domain [14]. The improved Pade approximation technique was combined with the modified pole clustering method in [15]. Various mixedmethod approaches were reported by employing the factor division technique along with the stability equation (SE) [16], modified pole clustering [17], eigen spectrum [18], and Routh approximation [19]. Finite-time reduction of linear systems using shifted Legendre polynomials was proposed in [20].
Optimal reduction of LTI SISO systems using the Luus-Jakola algorithm was reported in [21]. Metaheuristics such as the Cuckoo Search Algorithm (CSA) [22], modified CSA [23], Ant Lion Optimization (ALO) [24], and accelerated Opposition-based ALO (OB-ac-ALO) [24] algorithms were used to optimally determine the coefficients of the ROM by minimizing the integral square error (ISE) between the original system and the ROM, when subjected to a unit step input.
In [25]- [27], Big Bang Big Crunch (BBBC) optimization algorithm was used to determine the numerator polynomial of the ROM. On the other hand, the denominator polynomial was obtained using the SE [25], Routh approximation [26], and time moment matching method [27]. The denominator polynomial of the ROM was obtained by the SE technique, whereas Particle Swarm Optimization (PSO) and CSA were used to determine the optimal numerator coefficients in [28] and [29], respectively. In [30], the numerator and denominator polynomials of the ROM were obtained by using CSA and modified clustering technique, respectively. In [31], a generalized pole-clustering method determined the denominator coefficients of the ROM, whereas the numerator polynomial was determined using the mathematical approach presented in [8]. The numerator and denominator polynomials of the ROM were determined using the Jaya optimization algorithm and the eigen permutation method, respectively, in [32].
In [33], several nature-inspired computation techniques such as Firefly Algorithm (FA), GA, PSO, Differential Evolution (DE), Artificial Bee Colony, Harmony Search Algorithm (HSA), Invasive Weed Optimization, Grey Wolf Optimization, and Hybrid Grey Wolf Optimization (HGWO) algorithms were employed to determine the ROM based on pseudo-random binary sequence response matching (PRB-SRM). In [34], a multi-objective optimization framework based on the Routh-Pade approximation and HSA was reported for MOR design. Applications of the GA to solve MOR problems for single and multivariable systems were demonstrated in [35] and [36], respectively. ROM based on minimization of the ISE using ALO algorithm was also reported in [37].
The two main limitations of the existing research are outlined below: (i) The techniques reported in [22]- [24], and [37] determine the optimal coefficients of the ROM directly by minimizing the ISE based on the unit step response, which is a time-domain performance metric. Since steady-state matching condition is not incorporated, these methods may exhibit finite steady-state error (for e.g., [24], [37]). Furthermore, as reported in [33], the fast transient characteristics of the original system may not be clearly identified by employing the step input for MOR problems. (ii) The approaches reported in [25]- [32] fall under the category of 'mixed methods' which comprises two steps. In these techniques, the numerator and the denominator terms of the ROM are separately determined using two distinct methods. The numerator polynomial may be optimally approximated using the metaheuristics, whereas a non-optimal numerical technique is employed to determine the denominator coefficients. The timedomain characteristics of the original system are also considered here for the determination of the ROM.
Certain control system design tasks, such as loop shaping control [38], [39] may require an improved model reduction performance according to frequency-domain characteristics. Since the emphasis of the literature has primarily been on approximating the time-domain behavior of the full-order system, this paper presents an alternative approach to solve MOR problems. The contributions of this paper are the following: (i) A single-step, generalized, optimal, and frequencydomain-based MOR technique, which preserves the stability, minimum phase, and steady-state response, is presented. (ii) The suggested method is easily scalable, which can allow the designer to choose any desired order for the numerator and denominator polynomials of the ROM. (iii) The stability and minimum phase characteristics can be attained by setting a positive lower bound value for the design variables. This is simpler than [34], where inequality constraints based on Routh's criteria have to be satisfied to yield a stable model. (iv) The proposed approach can optimally obtain both the numerator and denominator coefficients of the ROM, which is in contrast to the reported mixed methods [25]- [32], where the optimizer determines the numerator dynamics only. Thus, the proposed technique is more straightforward and a real-parameter, single-objective, unconstrained optimization algorithm can implement the design procedure.
Odd Even k(s+z 0 ) Error (APE) are considered. (vi) The robustness of the proposed technique is demonstrated for all the examples through the convergence behavior of the algorithm and statistical analysis with respect to the quality of the solution and the ROM coefficients. The rest of the paper is organized as follows. In Section II, the proposed method is presented. Results of the MATLABbased experiments are discussed in Section III. The main conclusions of this work are outlined in Section IV.

II. PROPOSED TECHNIQUE
Without loss of generality, the transfer function of the ROM, as given by (2), can be redefined and classified into four categories based on the designers' requirement regarding the choice of M and N. Table 1 presents the proposed transfer function representation of G R (s) as a product of first-order and second-order functions depending on whether M and N take even or odd values. It may be noted that k, z i ,z i , p i , and p i (i = 0, 1, . . .) are the scalar parameters of the ROM. For the special case when Let the frequency response of the original (full-order) system and the proposed ROM be denoted by G O (jΩ) and G R (jΩ),respectively, where, Ω is the angular frequency expressed in radians/second (rad/s). The MOR design problem can be formulated in the optimal sense by minimizing the error in the frequency response between the full-order system and the ROM so that G R (jΩ) ≈ G O (jΩ).
Two essential requirements of a ROM for control system applications, such as design stability and minimum phase response, can be achieved by restricting the poles and zeros, respectively, of the model to lie on the left-half of the s-plane. While design stability guarantees the physical realization, the minimum phase behavior assures that the system exhibits the smallest group delay. For the ROMs proposed in Table 1, the stability and minimum phase characteristics can be achieved by restricting the design parameters (k, z i ,z i , p i ,p i ) to positive values only. Remark: To guarantee the stability of the ROM, the condition as defined by (3) needs to be satisfied: where i = 1, 2, . . ., N/2 (N : even), i = 0, 1, . . ., (N − 1)/2 (N : odd). (3) Proof: To attain a stable model in the s-domain, all the poles of the transfer function must lie in the left-half of the s-plane.
For N: even, the denominator polynomial can be written as: where Hence, the necessary and sufficient conditions to obtain a stable ROM are p i = (P i +P i ) andp i > 0. For N: odd, the additional factor (s+p 0 ) will always result in a left-hand side pole lying on the realaxis in the s-plane if p 0 > 0. A similar proof for the minimum phase response condition can be constructed utilizing the numerator polynomial of the proposed ROM transfer function. In an optimization procedure, the stability and minimum phase requirements can be met by setting the lower bounds of (k, z i ,z i , p i ,p i ) to be strictly greater than zero. Remark: For steady-state response matching, the condition as defined by (5) needs to be satisfied:

Steady-state matching condition
Proof: The steady-state response of the ROM approximates to steady-state of the original system when the following condition is satisfied.
By using the final value theorem, which suggests lim t→∞ f (t) = lim s→0 sF (s), one can write the following: This equation can be reorganized as lim s→0 (sG R (s) = sG 0 (s)) and written as: , the steady-state approximation condition as defined by (5) can be obtained.
Thus, the steady-state response of the original system for the four design cases can be achieved by the proposed ROM as per the conditions presented in Table 1. This condition is based on definingp N/2 (for N: even) andp (N −1)/2 (for N: odd) in terms of a 0 , b 0 ,p i , andz i . Such a condition can also be satisfied by the optimizer without employing any equality constraints. The proposed steady-state response matching condition results in the total number of design variables being one less than the total number of coefficients of the ROM.
The objective function for the proposed unconstrained optimization problem is defined by (9): where P is the total number of sampled frequency points that are logarithmically spaced in the interval Ω ∈ [Ω min , Ω max ] rad/s, and X denotes the vector of decision variables that needs to be determined by the optimizer. The vector X is defined in Table 1 for the four cases.
Since the utility of a ROM fundamentally lies in modeling the characteristics of the full-order system by employing small values of M and N, the reported works have primarily considered N = 2. Consequently, two design cases for the proposed model arise, such as {M = 0, N = 2} and {M = 1, N = 2}. The condition for steady-state response matching, the transfer function of the ROM, and the corresponding vector of decision variables (X) based on the presented approach for N = 2 are shown in Table 2.
It can be observed from Table 2 that both cases result in the same steady-state response matching criteria. Since this condition can be directly substituted in G R (s), consequently (i) the need to satisfy the equality condition using a constrained optimization technique is avoided, and (ii) the total number of decision variables is one less than the total number of coefficients of G R (s).
Based on the proofs presented above, the proposed optimization problem can be solved with guaranteed stability, minimum phase response, and steady-state response matching using a single-objective, real-parameter, unconstrained optimization algorithm by setting a positive value for the lower bound of X. In contrast to the mixed methods, the proposed technique optimally and simultaneously determines both the numerator and the denominator coefficients of the ROM transfer function.

III. SIMULATION EXPERIMENTS AND DISCUSSIONS
To quantify the accuracy of the models, various frequency and time response performance metrics as defined by (10)- (19) are employed. The smaller the value of the error terms as given by (10)-(18), the more accurate is the modeling accuracy. For exact matching of the steady-state response, the expression as defined by (19) should be equal to (a 0 /b 0 ) for the ROM.
Mean AM E = 20log 10 1 where y O (t) and y R (t) represent the step responses of the full-order system and the ROM, respectively, and L denotes the total number of logarithmic spaced sample points. An additional metric, namely the squared H 2 -norm H 2 2 of the difference between the transfer functions of the original and reduced systems, is also employed. It may be recalled that H 2 2 represents the ISE corresponding to an impulse input [40].
A state-of-the-art DE algorithm, namely the Enhanced Fitness-Adaptive Differential Evolution (EFADE) [41], is employed here as an optimization tool. EFADE employs a triangular mutation strategy to enhance the diversification and intensification characteristics as well as improve the convergence rate. The mutation operator consists of the convex combination vector of triplets based on three randomly chosen vectors and the difference vectors between the best, better, and worst among the selected vectors. In addition, the DE/rand/1/bin strategy is also utilized in EFADE. The adaptations of the scaling factor (F) and the crossover rate are carried out without introducing any new control parameters. In this work, the value of F is set to 0.7, which is in accordance with the specification reported in [41]. The minimization of the single-objective function defined by (9) using EFADE is carried out using a population size of 100, while the termination condition of the algorithm is set as 10000 numbers of function evaluations (NFEs).  . (20) The transfer function and pole-zero locations of the proposed ROM and the reported ones [9], [32], [34] are presented in Table 3. The zero {-5.0682} and poles {-0.0952, -50.1083} of the proposed model lie on the real-axis in the left-half of the s-plane, thus, exhibiting a stable and minimum phase behavior. In contrast, it may be noted that [9] yields a zero at 1.3301, which implies a non-minimum phase response, although the full-order system is of minimum phase type.
Detailed comparisons about the modeling accuracy of the proposed technique with the published literature are shown in Table 4. The results highlight that the proposed ROM outperforms [9], [32], [34] with respect to all error metrics. The ROMs reported in [9], [34] exhibit a finite steady-state error, while the proposed method alleviates this issue. Comparisons of magnitude and phase responses for the proposed model with the most competitive design reported in the literature (i.e., [32]) are illustrated in Figure 1(a). The reported model suffers a large deviation from the original systems' response regarding the magnitude and phase beyond approximately 1 rad/s. In contrast, the proposed model achieves significantly better response in the low as well as the high frequency range, as justified by the lower mean and max AME and APE metrics. This may be attributed to the placement of the zero {-5.0682} and one of the poles {-50.1083} of the proposed model being significantly different from that of the cited literature. Figure 1(b) shows the time-domain responses of the proposed ROM as compared to [32].    9.94 × 10 −4 } as compared to {0.06056, 6.73 × 10 −4 , 0.3045, 2.53 × 10 −3 } yielded for [32]. Therefore, the proposed method outperforms the recent literature on both the timeand frequency-domain performance indices.
Comparisons with the literature about the performance indices are presented in Table 6. The proposed ROM achieves the least values for all the frequency-domain error metrics. In particular, a significant improvement is attained regarding the mean AME for the proposed model (-26.67 dB), compared to the best-performing reported model [27] (-17.14 dB). It may be observed from Figure 2(a) that the proposed ROM attains better proximity to the original systems' frequency characteristics as compared to the cited literature, particularly in the interval from (10.14, 1000) rad/s and (4.24, 150) rad/s for magnitude and phase, respectively. Figure 2(b) presents the comparison of the time-domain responses with [27], which reveals that: (i) the impulse response of the proposed ROM achieves better agreement with that of the full-order systems' behavior. This superiority may also be confirmed from the comparison of H 2 2 metric between the proposed (0.0109) and the reported ROM [27] (0.3193), and (ii) the step response of the EFADE-based ROM and [27] is similar to the large-scale model. Quantitatively, the proposed ROM attains an improved performance over [27] about the ISE (0.001158 vs. 0.003916) and ITSE (0.004172 vs. 0.005397). Although [33] attains better H 2 2 (0.0356) as compared to [27] (0.3193), however, the step response performance measures of [33] are significantly inferior to [27] (see Table 6). Hence, the time responses comparison with [27] has been only considered in Figure 2(b). Overall, the proposed MOR technique achieves an improved modeling performance relative to the full-order system compared to the recently published literature.
G O (s) = 28s 3 + 496s 2 + 1800s + 2400 2s 4 + 36s 3 + 204s 2 + 360s + 240 . Table 7 presents the proposed and reported ROMs along with their pole-zero locations. It can be seen that: (i) the proposed model has real poles {-1.2121, -8.1797}, whereas all the reported ones provide complex-conjugate poles; and (ii) the zero for the proposed ROM is located at s = −6.5830, which implies that the poles and zeros are interlaced on the negative real-axis in the s-plane. The significant distinction in the pole-zero locations of the proposed model with those of [15], [30], [32] can be attributed to the fact that the cited literature aims to approximate the time-domain characteristics of the original system. Table 8 shows the comparison with the reported designs regarding the performance metrics. It is revealed that: (i) the proposed model outperforms all the reported ROMs about the various frequency-domain error indices. In particular, a sig-VOLUME 9, 2021    (ii) the IAE and ISE indices achieved by the proposed ROM are 0.6907 and 0.1110, respectively, which is superior to the best performances reported in the literature (IAE = 0.7179 [30] and ISE = 1.702 [15]); (iii) although the proposed approximant yields an inferior ITAE and ITSE, however, the H 2 2 achieved by the proposed 0.261 is significantly better than [15] 3.217, [30] 2.769, and [32] 3.819; and (iv) a finite steady-state error is exhibited by [30] since SSV = 9.999, while the proposed ROM attains the same SSV (= 10) as that of the original system.
As illustrated in Figure 3(a) (top), the magnitude plot of the proposed ROM is in good agreement with the original sys-tems' response throughout the design bandwidth. In contrast, the magnitude response of [30] is unable to approximate the ideal response beyond 4 rad/s accurately. As shown in Figure 3(a) (bottom), the phase plot of [30] also exhibits poor accuracy in the frequency interval of (1.4, 33.7) rad/s. In contrast, the phase response of the proposed ROM achieves better proximity to the full-order systems' behavior for the entire design range, as justified by a smaller max. APE (5.04 • ). The improved magnitude and phase responses of the proposed design may be possible due to the optimal interlacing of poles and zeros. From Figure 3(b), it can be observed that the impulse response of the proposed is distinctly superior when compared against [30], whereas the step response also attains a marginally better agreement with the original systems' behavior.

4) Example 4
The transfer function of the higher-order system given by (23) The transfer functions of the proposed ROM and the OBac-ALO algorithm-based model [24] are presented in Table 9. It may be emphasized that: (i) the proposed technique requires only two decision variables as compared to five needed by [24] to yield the ROM, and (ii) the steady-state response matching constraint is not incorporated in [24], which implies that the reported model may not exactly meet the SSV of the full-order system. Table 10 presents the comparisons with the cited literature for the various error metrics. Results reveal a marked improvement in both the time-and frequency-domain modeling characteristics of the proposed approximant. The improved accuracy of the EFADE-based model, when com-pared against [24], is also justified by the magnitude and phase-frequency responses, as shown in Figure 4(a). The magnitude and phase responses for [24] deviate from the original response beyond 1.933 rad/s and 0.292 rad/s, while the same for the proposed one are 13.78 rad/s and 2.719 rad/s. The superior time-domain responses of the proposed ROM can also be confirmed from Figure 4(b).

B. ROBUSTNESS ANALYSIS
The robustness of the ROMs designed using the proposed strategy for all the design examples is investigated in this section. Table 11 shows the minimum (min), max, mean, and standard deviation (SD) indices for the TE metric of the ROMs achieved using 30 independent trial runs of the proposed algorithm. The results highlight the extremely small (practically 0) value of SD about the TE, which implies that the same solution quality can be guaranteed regardless of the number of runs.   the ROM coefficients achieved using the proposed approach based on the same 30 runs of the algorithm. The design robustness is confirmed once again from the very small SD values attained for each of the model coefficients. Consequently, a single trial run of the algorithm is sufficient to generate the ROM using the proposed method. Further confirmation about the modeling robustness is demonstrated in Figure 5 by considering the convergence behavior of EFADE. It can be observed that the fitness convergence for a particular example for all the trial runs of the algorithm is the same. , where the worst-case NFE performance is 6100, 4500, 2900, and 8800, respectively. Hence, further reduction in t C , particularly for examples 1-3, is also possible to improve the design time.

C. GENERALITY OF THE PROPOSED APPROACH
In order to demonstrate the generic nature and scalability of the suggested strategy, the order diminution for the original system as defined by (21), i.e., Example 2, is considered here for approximation with third-order transfer function (N = 3). This particular example has been chosen since a     The transfer function of the third-order ROM reported in [15] is given by (25).
The error indices compared in Table 13 highlight the significant improvement of the proposed ROM in achieving superior frequency and time response characteristics as compared to [15]. Figures 6(a) and (b) show the frequency and time-domain response plots to confirm the improved modeling performance of the proposed technique graphically.

IV. CONCLUSIONS
An optimal, robust, and generalized framework for the order diminution of LTI SISO continuous-time systems is presented in this paper. The technique utilizes the frequencydomain information, unlike the published literature where the time-domain characteristic (e.g., ISE) of the full-order system is generally considered. The simplicity of the strategy as compared to the mixed methods is highlighted. The proposed model preserved the stability, minimum phase characteristic, and steady-state response of the original system. The improved accuracy over the recently published literature is extensively demonstrated using different performance indices. The excellent robustness and flexibility of the proposed scheme are also justified. Future work will extend the suggested scheme for reduction of fractional-order models [42] where the stability constraints may be implemented using the principal characteristic equations [43].