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

Received November 30, 2021, accepted December 11, 2021, date of publication December 13, 2021, date of current version December 23, 2021. Digital Object Identifier 10.1109/ACCESS.2021.3135279 A Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems SHIBENDU MAHATA 1, NORBERT HERENCSAR 2, (Senior Member, IEEE), BARIS BAYKANT ALAGOZ 3, AND CELALEDDIN YEROGLU 3 1Department of Electrical Engineering, Dr. B. C. Roy Engineering College, Durgapur, West Bengal 713206, India 2Department of Telecommunications, Faculty of Electrical Engineering and Communication, Brno University of Technology, 616 00 Brno, Czechia 3Department of Computer Engineering, Faculty of Engineering, Inonu University, 44000 Malatya, Turkey Corresponding author: Norbert Herencsar ([email protected]) This work was supported by the Project FEKT-S-20-6312. ABSTRACT 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. INDEX TERMS Frequency response, linear time-invariant system, model order reduction, reduced-order modelling, single-input-single-output system. I. INTRODUCTION Model order reduction (MOR) deals with the approximation of complex, high-order models with the low-order 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 generalized full-order (original) linear time-invariant (LTI) single-input-single-output (SISO) system is given by (1): GO(s)=amsm+am−1sm−1+ · · · + a1s+a0 sn+bn−1sn−1+ · · · + b1s+b0 ,(1) where ak(k=0,1,...,m) and bk(k=0,1,...,n−1) are the coefficients of the numerator and denominator, respectively, of GO(s); mand nare integers; and m≤n. The associate editor coordinating the review of this manuscript and approving it for publication was Zhiguang Feng . The transfer function of a ROM can be defined by (2): GR(s)=NR(s) DR(s)=AMsM+AM−1sM−1+ · · · + A0 sN+BN−1sN−1+ · · · + B0 ,(2) where Ak(k=0,1,...,M) and Bk(k=0,1,...,N−1) are the coefficients of the numerator and denominator polynomials, respectively, of GR(s); Mand Nare 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 VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ 165773 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems 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]. Geometric 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 LuusJakola 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 (PRBSRM). 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. (v) Four well-known design examples from the recently published literature are considered to exemplify the effectiveness of the proposed approach. Various performance measures pertaining to both the time domain, such as the Integral Absolute Error (IAE), ISE, Integral 165774 VOLUME 9, 2021 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 1. Transfer function of the proposed ROM for various combinations of Mand N. Time Absolute Error (ITAE), and Integral Time Square Error (ITSE), and the frequency domain, such as the maximum (max) and mean values about the Absolute Magnitude Error (AME) and Absolute Phase 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 MATLAB-based 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 Mand N. Table 1presents the proposed transfer function representation of GR(s) as a product of first-order and second-order functions depending on whether Mand N take even or odd values. It may be noted that k,zi,˜zi,pi, and ˜pi (i=0,1, . . .) are the scalar parameters of the ROM. For the special case when M=0,GR(s)=k/DR(s), where, k=A0. Let the frequency response of the original (full-order) system and the proposed ROM be denoted by GO(j) and GR(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 GR(j)≈GO(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,zi,˜zi,pi,˜pi) to positive values only. Remark: To guarantee the stability of the ROM, the condition as defined by (3) needs to be satisfied: pi=(Pi+˜ Pi),˜pi>0, 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: DR(s)= N/2 Y i=1 (s2+pis+ ˜pi)= N/2 Y i=1 (s+Pi)(s+˜ Pi),(4) where pi=(Pi+˜ Pi) and ˜pi=Pi˜ Pi. ...For stability, Pi>0 and ˜ Pi>0. Hence, the necessary and sufficient conditions to obtain a stable ROM are pi=(Pi+˜ Pi) and ˜pi>0. For N: odd, the additional factor (s+p0) will always result in a left-hand side pole lying on the real-axis in the s-plane if p0>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 VOLUME 9, 2021 165775 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 2. Transfer function of the proposed ROM (N=2). be met by setting the lower bounds of (k,zi,˜zi,pi,˜pi) to be strictly greater than zero. Remark: For steady-state response matching, the condition as defined by (5) needs to be satisfied: DR(s)|s=0=NR(s)|s=0 GO(s)|s=0 .(5) Proof: The steady-state response of the ROM approximates to steady-state of the original system when the following condition is satisfied. lim t→∞gR(t)=lim t→∞g0(t).(6) By using the final value theorem, which suggests lim t→∞f(t)= lim s→0sF(s), one can write the following: lim s→0sGR(s)=lim s→0sG0(s).(7) This equation can be reorganized as lim s→0(sGR(s)=sG0(s)) and written as: lim s→0(GR(s)=G0(s)) and GR(s)|s=0=G0(s)|s=0.(8) Then, substituting GR(s)=NR(s) DR(s), 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 defining ˜pN/2(for N: even) and ˜p(N−1)/2(for N: odd) in terms of a0,b0,˜pi, and ˜zi. 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): f=min (P X i=1 |GO(ji)−GR(ji,X)|),(9) where Pis the total number of sampled frequency points that are logarithmically spaced in the interval ∈ [min, max] rad/s, and Xdenotes the vector of decision variables that needs to be determined by the optimizer. The vector Xis defined in Table 1for 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 Mand 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 2that both cases result in the same steady-state response matching criteria. Since this condition can be directly substituted in GR(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 GR(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 (a0/b0) for the ROM. Mean AME =20log10 (1 L L X i=1 ||GO(ji)|−|GR(ji)||), [decibel (dB)] (10) Max AME =20log10 {max ||GO(j)|−|GR(j)||} ,[dB] (11) Mean APE =1 L L X i=1 |6GO(ji)−6GR(ji)|, [degrees (deg)] (12) Max APE =max |6GO(j)−6GR(j)|,[deg] (13) Total Error (TE)= L X i=1 |GO(ji)−GR(ji)|,(14) IAE =Z∞ 0 |yO(t)−yR(t)|dt,(15) 165776 VOLUME 9, 2021 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 3. Transfer function and pole-zero locations of the proposed and reported ROM for Example 1. ISE =Z∞ 0 [yO(t)−yR(t)]2dt,(16) ITAE =Z∞ 0 t|yO(t)−yR(t)|dt,(17) ITSE =Z∞ 0 t[yO(t)−yR(t)]2dt,(18) Steady State Value (SSV)=A0 B0 .(19) where yO(t) and yR(t) represent the step responses of the fullorder system and the ROM, respectively, and Ldenotes the total number of logarithmic spaced sample points. An additional metric, namely the squared H2-norm kH2k2of the difference between the transfer functions of the original and reduced systems, is also employed. It may be recalled that kH2k2represents 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 Fis 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). The lower and upper bounds for all decision variables are set as 10−6 and 200, respectively, P=100, and L=1000 (logarithm spacing) for all design examples. For n−m=1, the values of Mand Nare set to 1 and 2, respectively. On the other hand, for n−m>1, the values of Mand Nare chosen as 0 and 2, respectively. Thirty independent trial runs of the proposed optimization routine are carried out for each of the considered examples. All the simulations are conducted in the following environment – CPU: Intel i3 @ 1.70 GHz, RAM: 2 GB, Operating System: Windows 7 (64 bit), and Software: MATLAB 2014a. A. PERFORMANCE ANALYSIS 1) EXAMPLE 1 Consider a sixth-order system as defined by (20), where the zeros and poles are located at {–0.2632, –0.3530, –0.7415, –0.0711 ±2.6929i}and {–0.1, –0.2, –0.5, –1, –5, –10}, respectively. GO(s) =2s5+3s4+16s3+20s2+8s+1 2s6+33.6s5+155.94s4+209.46s3+102.42s2+18.3s+1. (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]. The impulse response of the designed ROM is closer to that of the original systems’ especially in the interval [0, 0.1] s. While the step responses for the proposed and the reported models are in agreement with the theoretical characteristics, the proposed ROM achieves superior accuracy as justified by {IAE,ISE, ITAE,ITSE}of {0.03824, 5.03 ×10−4, 0.1214, 9.94 ×10−4} as compared to {0.06056, 6.73 ×10−4, 0.3045, 2.53 ×10−3} VOLUME 9, 2021 165777 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 4. Comparison with the literature for Example 1. FIGURE 1. (a) Frequency and (b) time responses comparison plots with the reported literature [32] for Example 1. yielded for [32]. Therefore, the proposed method outperforms the recent literature on both the timeand frequency-domain performance indices. 2) EXAMPLE 2 Consider the eighth-order system as given by (21), as shown at the bottom of this page. The zeros and poles of the full-order system are located at {–0.3195, –5.8917, –7.4292, –2.4498 ±0.5278i, –5.0078 ±0.6543i}and {–1, –2, –3, –4, –5, –6, –7, –8}, respectively. Table 5shows the ROMs (proposed and reported) along with their corresponding pole-zero values. It may be noted that the ROM reported in [27] will yield a finite steady-state error since A0B0=1.0005 whereas the proposed and the original systems’ SSV(=1) are the same. 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 kH2k2metric 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 largescale 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 kH2k2(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 GO(s)=18s7+514s6+5982s5+36380s4+122664s3+222088s2+185760s+40320 s8+36s7+546s6+4536s5+22449s4+67284s3+118124s2+109584s+40320.(21) 165778 VOLUME 9, 2021 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 5. Transfer function and pole-zero locations of the proposed and reported ROM for Example 2. TABLE 6. Comparison with the literature for Example 2. FIGURE 2. (a) Frequency and (b) time responses comparison plots with the reported literature [27] for Example 2. an improved modeling performance relative to the full-order system compared to the recently published literature. 3) EXAMPLE 3 In this example, the transfer function of the original system as given by (22) has poles and zeros located at {–1.1967 ±0.6934i, –7.8033 ±1.3576i}and {–2.1612 ±1.3152i, –13.3919}, respectively. GO(s)=28s3+496s2+1800s+2400 2s4+36s3+204s2+360s+240.(22) Table 7presents 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 8shows the comparison with the reported designs regarding the performance metrics. It is revealed that: (i) the proposed model outperforms all the reported ROMs VOLUME 9, 2021 165779 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 7. Transfer function and pole-zero locations of the proposed and reported ROM for Example 3. TABLE 8. Comparison with the literature for Example 3. FIGURE 3. (a) Frequency and (b) time responses comparison plots with the reported literature [30] for Example 3. about the various frequency-domain error indices. In particular, a significant improvement is achieved in {mean AME, max APE}for the proposed {–24.55 dB, 5.04◦}as compared to the most competitive reported model [30] {–18.15 dB,20.34◦}; (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 kH2k2achieved 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 systems’ 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 165780 VOLUME 9, 2021 S. Mahata et al.: Robust Frequency-Domain-Based Order Reduction Scheme for Linear Time-Invariant Systems TABLE 9. Transfer function and pole-zero locations of the proposed and reported ROM for Example 4. FIGURE 4. (a) Frequency and (b) time responses comparison plots with the reported literature [24] for Example 4. 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) has a single zero at s= −4 and real, non-repeated poles located at s= { − 1,−3,−5,−10}. GO(s)=s+4 s4+19s3+113s2+245s+150.(23) The transfer functions of the proposed ROM and the OB-ac-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 steadystate 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 timeand frequency-domain modeling characteristics of the proposed approximant. The improved accuracy of the EFADE-based model, when compared 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. Table 12 shows the statistical indices for the values of 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 5by 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. The {min, max, mean, SD} values of computational time tC(expressed in seconds) for examples 1, 2, 3, and 4 with 10000 NFEs are {2.165, 2.627, 2.207, 0.0832},{2.837, 3.266, 2.920, 0.0916},{2.298, 2.399, 2.348, 0.0197}, and {1.522, 1.734, 1.549, 0.0349}, VOLUME 9, 2021 165781