Iterative Rational Krylov Algorithms for model reduction of a class of constrained structural dynamic system with Engineering applications
Xin Du, M. Monir Uddiny, A. Mostakim Fonyz, Md. Tanzim Hossainx, Md. Nazmul Islam Shuzan
aa r X i v : . [ m a t h . O C ] J a n Iterative Rational Krylov Algorithms for modelreduction of a class of constrained structuraldynamic system with Engineering applications
Xin Du ∗ , M. Monir Uddin † , A. Mostakim Fony ‡ , Md. TanzimHossain § , and Md. Nazmul Islam Shuzan ¶ Abstract
This paper discusses model order reduction of large sparse second-order index-3 differential algebraic equations (DAEs) by applying Itera-tive Rational Krylov Algorithm (IRKA). In general, such DAEs arise inconstraint mechanics, multibody dynamics, mechatronics and many otherbranches of sciences and technologies. By deflecting the algebraic equa-tions the second-order index-3 system can be altered into an equivalentstandard second-order system. This can be done by projecting the systemonto the null space of the constraint matrix. However, creating the pro-jector is computationally expensive and it yields huge bottleneck duringthe implementation. This paper shows how to find a reduce order modelwithout projecting the system onto the null space of the constraint ma-trix explicitly. To show the efficiency of the theoretical works we applythem to several data of second-order index-3 models and experimentalresultants are discussed in the paper. keywords : Structured index-3 differential algebraic equations, sparsity,Model order reduction, Iterative Rational Krylov Algorithms. ∗ School of Mechatronic Engineering and Automation, Shanghai University, Shanghai-200072, China and Key Laboratory of Modern Power System Simulation and Control &Renewable Energy Technology, Ministry of Education(Northeast Electric Power University),Jilin-132012, China , [email protected] † Department of Mathematics and Physics, North south University, Dhaka-1229,Bangladesh, [email protected] ‡ Department of Mathematics, Chittagong University, Chittagong, Bangladesh, [email protected] § Department of Electrical and Computer Engineering, North South University, Dhaka-1229, Bangladesh, [email protected] ¶ Department of Electrical and Computer Engineering, North South University, Dhaka-1229, Bangladesh, [email protected] Introduction
In mechanics or multibody dynamics linearized equation of motion with holo-nomically constraint has the following form [1, 2] M ¨ x ( t ) + D ˙ x ( t ) + Kx ( t ) + G T z ( t ) = Fu ( t ) , (1a) Gx ( t ) = 0 , (1b) y ( t ) = Lx ( t ) , (1c)where M , K , D ∈ R n × n are sparse matrices known as mass, stiffness anddamping matrices respectively. x ( t ) ∈ R n , u ( t ) ∈ R m and y ( t ) ∈ R q arerespectively known as states, inputs and outputs vectors. The constraint matrix G ∈ R n × n (with n < n ) is associated with the given algebraic constraints z ( t ) ∈ R n . Furthermore, F ∈ R n × m is the input matrix corresponding tothe input vector u ( t ) and L ∈ R q × n is the output matrix associated to themeasurement output vector y ( t ). Such structured dynamical system also appearin mechatronics where electrical and mechanical parts are coupled or in theelectric circuits [3].If we convert the system into first-order form then it becomes first-orderindex-3 system [2]. Therefore the system in (1) is called second-order index-3descriptor system. If the system becomes very large then it is very expensiveto simulate, control and optimize. Therefore, we want to reduce the complexityof the model through model order reduction (MOR) [4, 5, 6]. Among differentMOR methods [4, 5, 6] the two most frequently applied modern MOR meth-ods are the balanced truncation (BT) [7] and the rational interpolation of thetransfer function by the iterative rational Krylov algorithm (IRKA) [8]. Bothapproaches have been extended to first-order descriptor systems [9, 10]. The bal-ancing based model order reduction of second-order index-3 system (1) has beeninvestigated for second-order-to-first-order and second-order-to-second-order re-ductions in [2] and [11] respectively. On the other hand, the authors in [12]discussed IRKA for the model reduction of the underlying descriptor system.In order to follow the proposed algorithm one has to convert the system into afirst-order form. Besides at each iteration one has to solve a linear system withdimension 2( n + n ) which is computationally expensive tasks.In this paper we discuss second-order-to-second-order model reduction ofsecond-order index-3 descriptor system via IRKA without converting the systeminto first-order form. In the literature second-order-to-second-order reductionis called structure preserving model order reduction (SPMOR). IRKA basedSPMOR for the standard second-order system was developed by Wyatt in hisP.hD., thesis [13]. This idea was generalized for the second-order index-1 systemwhich is slightly different from (1) in [14]. Like index-1 system the second-orderindex-3 system can be converted into standard second-order system. In thiscase instead of using Schur complement techniques as used in [14] we applyprojection onto hidden manifold. This idea was already found in [10, 9, 15] forthe firs-order index-2 systems. On the other hand, for the second-order index-3 system the technique was implemented using a balancing based model orderreduction. However, there was no investigation of this idea for this system usingthe IRKA. This paper contributes to close this gap. That is we mainly devoteto second-order-to-second-order model reduction of second-order index-3 systemusing IRKA. Following the procedure in [10, 9] first we show that the second-2rder index-3 descriptor system can be projected onto the null space of theconstraint matrix which we call hidden manifold to obtain a standard second-order system. Then we can apply the technique as in [13] to obtain a standardsecond-order reduced order system. It is shown in the paper that the explicitcomputation of hidden manifold projector is not required. This is importantbecause creating a hidden manifold projector demands a lot of computationaltimes. Moreover, the projected system is converted into a dense form whichyields huge bottleneck in implementing the reduced order model. The proposedmethod is applied to several models coming from Engineering applications. Theperformances of the proposed algorithm seems to be promising and the resultsare better than balancing based techniques in both approximation accuracy andcomputational time which appears in numerical results.Rest of the article is organized as follows. Section 2 briefly discuss IRKAbased SPMOR of second-order system and reformulation of second-order index-3 system from previous literature which are the main ingredient to obtain thenew results of this paper. The main contribution of this paper will be discussedin Section 3. In this section we developed IRKA based SPMOR of second-orderindex-3 system. The subsequent section illustrates numerical results. At theend, Section 5 presents the conclusive remarks. In the following texts at first we briefly discuss IRKA for standard second-ordersystem to obtain second-order reduced order model. Then we will show howto convert the second-order index-3 system into second-order standard systemby projecting onto the null space of the constraint matrix. In fact we establishsome definitions and notations based on the previous literature that will be inthe upcoming sections.
Structure preserving IRKA (SPIRKA) for a second-order standard system wasproposed in [13]. The SPIRKA is mainly based on the IRKA of first-ordersystem which was originally proposed in [8]. This prominent algorithm wasdeveloped by Gugercin et al., in [8] to achieve the H -optimal model reductionvia interpolatoy projection technique. To explain the SPIRKA let us considera second-order linear time-invariant (LTI) continuous-time system M ¨ ξ ( t ) + D ˙ ξ ( t ) + K ξ ( t ) = F u ( t ) ,y ( t ) = L ξ ( t ) , (2)where M , D and K are non-singular, and ξ ( t ) is the n dimensional state vector.Consider that the system is MIMO and its transfer function is defined by T ( s ) = L ( s M + s D + K ) − F ; s ∈ C . (3)Our goal is to obtain an r dimensional ( r ≪ n ) reduce order modelˆ M ¨ˆ ξ ( t ) + ˆ D ˙ˆ ξ ( t ) + ˆ K ˆ ξ ( t ) = ˆ F u ( t ) , ˆ y ( t ) = ˆ L ˆ ξ ( t ) , (4)3here the reduced coefficient matrices are constructed asˆ M = W T MV , ˆ D = W T DV , ˆ K = W T KV , ˆ F = W T F , ˆ L = LV , (5)and the transfer function of the reduced order model can be defined asˆ T ( s ) = ˆ L ( s ˆ M + s ˆ D + ˆ K ) − ˆ F ; s ∈ C . (6)According to [13] the procedure of IRKA for second-order system is same asthe first-order system. We want to construct reduced order model (4) in suchway that the reduced transfer function (6) interpolate to the original transferfunction (3) at some interpolation points. Moreover the reduced order modelsatisfies the interpolation conditions mentioned below [13].Given a set of interpolation points { α , α , · · · , α r } ⊂ C , and sets of leftand right tangential directions { c , c , · · · , c r } ⊂ C m , { b , b , · · · , b r } ⊂ C p arerespectively defined by V = [ v F b , v F b , · · · , v r F b r ] ,W = (cid:2) w L T c , w L T c · · · , w r L T c r (cid:3) , (7)Where v i = ( α i M + α i D + K ) − and w i = ( α i M T + α i D T + K T ) − ; i = 1 , , · · · r .If the reduced-order model (4) is constructed by V and W , the reduced transfer-function (6) tangentially interpolates (3), satisfies the interpolation conditions T ( α i ) b i = ˆ T ( α i ) b i ,c Ti T ( α i ) b i = c Ti ˆ T ( α i ) b i ,c Ti T ′ ( α i ) b i = c Ti ˆ T ′ ( α i ) b i , (8)for i = 1 , , . . . , r, which is known as Hermite bi-tangential interpolation con-ditions. One of the challenging parts of SPIRKA is to find a set of optimalinterpolation points as well as tangential directions since they are not prede-fined. In [13] author shows several remedies. Among them this paper considerthe following strategy. Construct ( E, A, B, C ) := (cid:18)(cid:20) I
00 ˆ M (cid:21) , (cid:20) I − ˆ K − ˆ D (cid:21) , (cid:20) F (cid:21) , (cid:2) ˆ L (cid:3)(cid:19) . (9) Then find r dimensional reduced order model ( ˆ E, ˆ A, ˆ B, ˆ C ) from ( E, A, B, C ).The interpolation points and tangential directions for the next iteration stepare constructed from the mirror image of the eigenvalues and the eigenvectorsof (
A, E ). The reduced order model ( ˆ E, ˆ A, ˆ B, ˆ C ) can be constructed again bythe IRKA. Note that IRKA of first-order system is presented in [6, Algorithm1]. The whole procedure for SPIRKA is summarized in Algorithm 1. We already have mentioned in earlier section that projecting the index-3 system(1) onto the hidden manifold we can convert the system into an index-0 i.e.,second-order standard system like (2). However, such conversion for a large4 lgorithm 1:
SPIRKA for standard Second-Order Systems.
Input : M , D , K , H , L . Output: ˆ M , ˆ D , ˆ K , ˆ H , ˆ L . Consider: interpolation points { α i } ri =1 ⊂ C , left and right tangentialdirections { b i } ri =1 ⊂ C p and { c i } ri =1 ⊂ C m . Form V = [ v F b , v F b , · · · , v r F b r ] ,W = (cid:2) w L T c , w L T c · · · , w r L T c r (cid:3) , Where v i = ( α i M + α i D + K ) − and w i = ( α i M T + α i D T + K T ) − ; i = 1 , , · · · r . while (not converged) do ˆ M = W T MV , ˆ D = W T DV , ˆ K = W T KV , ˆ F = W T F , ˆ L = LV . Use the first-order representation (
E, A, B, C ) as in (9) find thereduced-order matrices ˆ E , ˆ A , ˆ B and ˆ C . Compute ˆ Az i = λ i ˆ Ez i and y ∗ i ˆ A = λ i y ∗ i ˆ E for α i ← − λ i , b ∗ i ← − y ∗ i ˆ B and c ∗ i ← ˆ Cz ∗ i for all i = 1 , · · · , r Repeat Step 2. i = i + 1; end while Construct the reduced matrices by repeating Step 4scale dynamical system is practically impossible due to additional complexities.The idea of conversion is already developed in [11]. For our convenience, webriefly introduce this in the following.Let us consider the projector onto the null-space of G ,Π := I − G T ( GM − G T ) − GM − , (10)which satisfies Π M = M Π T , Null(Π) = Rang( G T ), Rang(Π) = Null( GM − )and the most importantly Gx = 0 ⇐⇒ Π T x = x . (11)Readers are referred to e.g., [11] to see the details of these properties with proofs.Now applying these identities into (2) we obtainΠ M Π T ¨ x ( t ) + Π D Π T ˙ x ( t ) + Π K Π T x ( t ) = Π F u ( t ) , (12a) y ( t ) = L Π T x ( t ) . (12b)The dynamical system (12) still has unnecessary equations due to the singularityof Π. Those equations can be avoided by splitting Π = Ψ l Ψ r , where Ψ l , Ψ r ∈ R n × ( n − n ) and they satisfies Ψ Tl Ψ r = I n − n , (13)where I n − n is an ( n − n ) × ( n − n ) identity matrix. Inserting the decom-position of Π into (12) and considering ˜ x ( t ) = Ψ lT x ( t ), the resulting dynamical5ystem leads toΨ Tr M Ψ r ¨˜ x ( t ) + Ψ Tr D Ψ r ˙˜ x ( t ) + Ψ Tr K Ψ r ˜ x ( t ) = Ψ Tr F u ( t ) , (14a) y ( t ) = L Ψ r ˜ x ( t ) . (14b)This system is now a standard second-order system as described by (2). In factsystem in (14) can be seen as the system (12) with the redundant equationsbeing removed through the Ψ r projection. Note that the coefficient matrices of(14) are dense if compared to (1). Therefore, for the large-scale index-3 systemexplicit computation of (14) is forbidden. In fact the dynamical systems (1),(12) and (14) are equivalent in the sense that they are the different realizationsof the same transfer function. Moreover, their finite spectra are the same whichhas been proven in the sequel. Once the index-3 system (1) is converted intothe index-0 system (14) then the SPIRKA i.e., Algorithm 1 can be applied tothe converted system. However, as the converted system is dense, the computa-tional costs of Algorithm 1 becomes high. For a very high dimensional index-3descriptor system computing (14) is not possible due to the restriction of com-puter memory. Therefore we are motivated to construct reduced order modelwithout forming system (14) explicitly. The SPIRKA introduced in Section 2 can be applied to the projected sys-tem (14). As already mentioned, this is infeasible for a large-scale system.Therefore, the technique can be applied to the equivalent system (12) instead.For this purpose, following the discussion in Section 2, we can create the rightand left projectors as V = [ v , v , · · · , v r ] ,W = [ w , w · · · , w r ] , (15)where v i = ( α i ˜ M + α i ˜ D + ˜ K ) − ˜ F b i ,w i = ( α ˜ M T + α ˜ D T + ˜ K T ) − ˜ L T c i , for i = 1 , , · · · r and in which ˜ M = Π M Π T , ˜ D = Π D Π T , ˜ K = Π K Π T , ˜ F =Π F Π T and ˜ L = L Π T . The main expensive task here is to compute each vectorinside the projectors by solving a linear system. For example to construct V ,at i -th iteration we find ( α i ˜ M + α i ˜ L + ˜ K ) − ˜ F b i to solve the linear system( α i ˜ M + α i ˜ D + ˜ K ) v i = ˜ F b i , which implies Π( α i M + α i D + K )Π T v i = F b i . (16)This linear system can be solved efficiently by applying the following Lemma .6 emma 3.1. The matrix κ satisfies κ = Π T κ and Π( α M + α D + K )Π T κ =Π F b if and only if (cid:20) α M + α D + K G T G (cid:21) (cid:20) κ Λ (cid:21) = (cid:20) F b (cid:21) . (17) Proof. If κ = Π T κ , then by using (11) we have G κ = 0 (18)which is the second block of equation (20). Furthermore Π( α M + α D + K )Π T κ = Π F b implies Π( α M + α D + K ) κ = Π F b, Π(( α M + α D + K ) κ − F b ) = 0 . This means ( α M + α D + K ) κ − F b is in the null space of Π. We know thatNull(Π) = Rang( G T ). Therefore, there exists Λ such that ( α M + α D + K ) κ − F b = − G T Λ which implies( α M + α D + K ) κ + G T Λ = F b. (19)Equations (18) and (19) yield (20). Conversely, we assume (20) holds. Fromthe second line of (20) we obtain G κ = 0 and thus κ = Π T κ . Now from firstequation we obtain ( α M + α D + K ) κ + G T Λ = F b. Applying (11) this equation gives( α M + α D + K )Π T κ + G T Λ = F b. Multiplying both sides by Π and since Π G T = 0 we haveΠ( α M + α D + K )Π T κ = Π F b. This completes the proof.Following Lemma 3.1 instead of solving (16) we can solve (cid:20) α i M + α i D + K G T G (cid:21) (cid:20) v i Λ (cid:21) = (cid:20) F b i (cid:21) , (20)for ξ . Although this system is larger than its projected system. Thereforewe solve this linear system to avoid constructing the projector. Similarly toconstruct W at each iteration we compute w i = ( α i ˜ M T + α i ˜ D T + ˜ K T ) − ˜ L T c i by solving the following linear system (cid:20) α i M T + α i D T + K T G T G (cid:21) (cid:20) w i Γ (cid:21) = (cid:20) L T c i (cid:21) . Once we have V and W , apply them to (12) to find the reduce order modelˆ M ¨ˆ x ( t ) + ˆ D ˙ˆ x ( t ) + ˆ K ˆ x ( t ) = ˆ F u ( t ) , (21a)ˆ y ( t ) = ˆ L ˆ x ( t ) , (21b)7n which the reduced matrices are constructed as followsˆ M = W T Π M Π T V, ˆ D = W T Π D Π T V, ˆ K = W T Π K Π T V, ˆ F = W T Π F , ˆ L = L Π T V. Due to the properties of the projector, as mentioned in (11) we have Π T V = V and Π T W = W or (Π T W ) T = W T . Therefore the reduced matrices can beconstructed without using Π as followsˆ M = W T M V, ˆ D = W T D V, ˆ K = W T K V, ˆ F = W T F , ˆ L = L V. The whole procedure to construct reduced order model from second-order index-3 system (1) is summarize in Algorithm 2.
Update interpolation points and tangential directions.
In IRKA weneed to update the interpolation points and tangential direction at each itera-tion steps which is often a challenging task. Usually, the interpolation pointsand tangential direction are updated by using mirror image of the eigenvaluesand eigenvector of the reduced order model. In SPIRKA this task is compli-cated because we need to solve the quadratic eigenvalue problems. If we solvea quadratic eigenvalue problem [16] using ( ˆ M , ˆ D , ˆ K ) we obtain 2 r number ofeigenvalues and eigenvectors. Selecting r number of optimal interpolation pointsand corresponding tangential direction is challenging task. To resolve this com-plexity [13] propose a techniques which is discussed in Algorithm 2. To followthis idea, we need to apply standard IRKA onto the converted first order sys-tem from the reduced second order system. However, this is again an iterativemethod which is computationally expensive. In this paper to construct theinterpolation points and tangential directions we construct E := (cid:20) I M (cid:21) , A := (cid:20) I − ˆ K − ˆ D (cid:21) ,B := (cid:20) F (cid:21) and C := (cid:2) ˆ L (cid:3) . (22)Then we apply MATLAB function balred to compute r dimensional reduced-order model from ( ˆ E, ˆ A, ˆ B, ˆ C ). The interpolation points are updated by choos-ing the mirror images of the eigenvalues of the pair ( ˆ A, ˆ E ) as the next interpo-lations points. This seems to be more efficient than the existing one. To asses the efficiency of the proposed algorithm, i.,e., Algorithm 2 we haveapplied this to two sets of data. First one is coming from a damped spring-masssystem (DSMS) with holonomic constraint which is taken from [17]. The secondset of data set is a constrain triple chain oscillator model (TCOM). This data isoriginated in [18] but with the index-3 setup described in [2]. The details of themodels is also available in [11]. We intentionally have considered the same datafor the both model examples as in [11] since we want to compare the results of8 lgorithm 2:
IRKA for Second-Order Index-3 Descriptor Systems.
Input : M , D , K , F , L . Output: ˆ M , ˆ D , ˆ K , ˆ F , ˆ L Select randomly a set of the interpolation points { α i } ri =1 and thetangential directions { b i } ri =1 and { c i } ri =1 . Construct the projection matrices V s = [ v , v , · · · , v r ] and W s = [ w , w , · · · , w r ] , where v i & w i ; i = 1 , · · · , r are the solutions of the linear systems (cid:20) α i M + α i D + K G T G (cid:21) (cid:20) v i Λ (cid:21) = (cid:20) F b i (cid:21) and (cid:20) α i M T + α i D T + K T G T G (cid:21) (cid:20) w i Γ (cid:21) = (cid:20) L T c i (cid:21) , respectively. while (not converged) do Constructˆ M = W T M V, ˆ D = W T D V, ˆ K = W T K V, ˆ F = W T F , ˆ L = L V. Compute (
E, A, B, C ) as in (22), then using MATLAB function balred compute r dimensions model ( ˆ E, ˆ A, ˆ B, ˆ C ). Compute ˆ Az i = λ i ˆ Ez i and y ∗ i ˆ A = λ i y ∗ i ˆ E for α i ← − λ i , b ∗ i ← − y ∗ i ˆ B and c ∗ i ← ˆ Cz ∗ i . Repeat Step 2. i = i + 1. end while Finally construct the reduced-order matricesˆ M = W T M V, ˆ D = W T D V, ˆ K = W T K V, ˆ F = W T F , ˆ L = L V. this paper with that one. The dimension of the models including the number ofdifferential and algebraic equations and inputs/outputs are displayed in Table 1.This experiment was carried out with MATLAB ® R2015a (8.5.0.197613) ona board with 4 × INTELCore TM i5-4460s CPU with a 2.90 GHz clock speed and16 GB RAM.We apply Algorithm 2 to the both models and find 30 dimensional reduced-order models. Figures 1 and 4 show the frequency domain analysis of the originaland reduced models for the DSMS and TCOM, respectively. The frequency re-sponses of the full and the reduced-order models and their absolute and relativeerrors for the DSMS are shown in Figure 1 over the frequency interval [10 − , ].Sub-figure 1a shows that the frequency responses of the reduced-order modelsare matching with the original model with good accuracy. The absolute and rel-9odels dimension n and n inputs/outputsDSMS 2200 2000 and 200 1/3TCOM 11001 6001 and 5000 1/1Table 1: The dimension of the tested models including number of differentialand algebraic variables and, inputs and outputs.ative errors between full and the reduced-order models are shown in Sub-figures1b and 1c, respectively. On the other hand, Figure 4 depicts the frequencyresponses, absolute and relative errors of the full and the reduced-order modelsof the TCOM over the frequency interval [10 − , ]. This figure also shows(in Sub-figure 2a) that the frequency responses of the reduced-order modelsare matching correctly with the original model. Sub-figures 2b and 2c shows agood approximation between the original and reduced-order models using theabsolute and relative errors, respectively. Comparisons of Balanced truncation and IRKA.
To compare the per-formance of IRKA and balanced truncation we compute 30 dimensional reduced-order model applying [11, Algorithm 2] to the TCOM. This algorithm can com-pute several reduced-order models based on different balancing criterion. Herewe consider the velocity-velocity balancing label which gives the best approxima-tion. Figure 3 shows the approximation errors of 30 dimensional reduced-ordermodels computed by the BT and IRKA. From Figure 3 it seems that the per-formance of IRKA is better than the BT. Both the absolute error and relativeerrors as shown in Figures 3a and 3b, respectively, IRKA depicts better accu-racy than the balanced truncation. On the other hand, when we consider thecomputation time, again the performance of IRKA is far better than the BTwhich is reflected in Figure 4. We know that balanced truncation is expensivemethod since it requires to solve two continuous-time algebraic Lyapunov equa-tions. The solution of the Lyapunov equations involved the computation of shiftparameters which is a computational. We have solved the Lyapunov equationsby [11, Algorithm 3] using adaptive shift parameters. See, e.g., [11] for details.Note that the computational time of IRKA is increasing if the dimension ofreduced order model and the number of iterations are increased gradually.
In this paper we have discussed a IRKA based technique to find a reducedsecond-order system from a large-scale sparse second-order index-3 system. Inparticular, we have linearized equation of motion with holonomic constrainswhich arise in constrained mechanics or multibody dynamics. It has been shownthat the index-3 system can be converted into index-0 by projecting onto thehidden manifold to apply the standard second-order IRKA. But creating projec-tor is often computationally expensive task and it yields system matrices dense.Therefore we have modified the standard IRKA for the underlying index-3 de-scriptor system. We also have shown a clever techniques to compute the inter-polation points and tangential directions. The proposed algorithm was applied100 − − − − ω σ m a x ( T ( j ω )) full reduced model (a) Frequency response. − − − − − ω σ m a x ( T ( j ω ) - ˆ T ( j ω )) (b) Absolute error. − − − − − ω σ m a x ( T ( j ω ) − ˆ T ( j ω )) σ m a x ( T ( j ω )) (c) Relative error. Figure 1: Comparison of original and the 30 dimensional reduced models forthe DSMS.to several data of second-order index-3 models. Numerical results showed thatthe proposed algorithm can generated lower dimensional model with higher ac-curacy. The IRKA based method is better than Balanced truncation in termsof accuracy and computational complexity as well.
This research work was funded by NSU-CTRG research grant under the projectNo.: CTRG-19/SEPS/05. It was also supported by National Natural ScienceFoundation of China under Grant No. (61873336, 61873335), the FundamentalResearch Funds for the Central Universities under Grant (FRF-BD-19-002A),and the High-end foreign expert program of Shanghai University,110 − − − ω σ m a x ( T ( j ω )) full reduced model (a) Frequency response. − − − − − − ω σ m a x ( T ( j ω ) - ˆ T ( j ω )) (b) Absolute error. − − − − − − ω σ m a x ( T ( j ω ) − ˆ T ( j ω )) σ m a x ( T ( j ω )) (c) Relative error. Figure 2: Comparison of the original and 30 dimensional reduced models forthe TCOM. 120 − − − − − − ω σ m a x ( T ( j ω ) - ˆ T ( j ω )) IRKA BT (a) Absolute error. − − − − − − ω σ m a x ( T ( j ω ) − ˆ T ( j ω )) σ m a x ( T ( j ω )) (b) Relative error. Figure 3: Comparison of the original and 30 dimensional reduced models com-puted by IRKA and balanced truncation for the TCOM.
81% 19%
BTIRKA
Figure 4: Time comparisons of both balanced truncation and IRKA for theTCOM. 13 eferences [1] E. Eich-Soellner and C. F¨uhrer,
Numerical Methods in Multibody Dynamics ,ser. European Consortium for Mathematics in Industry. Stuttgart: B. G.Teubner GmbH, 1998.[2] M. M. Uddin, “Gramian-based model-order reduction of constrained struc-tural dynamic systems,”
IET Control Theory & Applications , vol. 12, no. 7,p. 2337 – 2346, 2018.[3] R. Riaza,
Differential-Algebraic Systems. Analytical Aspects and CircuitApplications . Singapore: World Scientific Publishing Co. Pte. Ltd., 2008.[4] A. Antoulas,
Approximation of Large-Scale Dynamical Systems , ser. Ad-vances in Design and Control. Philadelphia, PA: SIAM Publications,2005, vol. 6.[5] F. Bennini, “Ordnungsreduktion von elektrostatisch-mechanischen FiniteElemente Modellen auf der Basis der modalen Zerlegung,” Ph.D. Thesis,Technische Universit¨at Chemnitz, Chemnitz, 2005.[6] M. M. Uddin,
Computational Methods for Approximation of Large-ScaleDynamical Systems . New York, USA: Chapman and Hall/CRC, 2019.[7] B. C. Moore, “Principal component analysis in linear systems: controlla-bility, observability, and model reduction,”
IEEE Trans. Autom. Control ,vol. AC–26, no. 1, pp. 17–32, 1981.[8] S. Gugercin, A. C. Antoulas, and C. A. Beattie, “ H model reduction forlarge-scale dynamical systems,” SIAM J. Matrix Anal. Appl. , vol. 30, no. 2,pp. 609–638, 2008.[9] M. Heinkenschloss, D. C. Sorensen, and K. Sun, “Balanced truncationmodel reduction for a class of descriptor systems with application to theOseen equations,”
SIAM J. Sci. Comput. , vol. 30, no. 2, pp. 1038–1063,2008.[10] S. Gugercin, T. Stykel, and S. Wyatt, “Model reduction of descriptor sys-tems by interpolatory projection methods,”
SIAM J. Sci. Comput. , vol. 35,no. 5, pp. B1010–B1033, 2013.[11] M. M. Uddin, “Structure preserving model order reduction of a class ofsecond-order descriptor systems via balanced truncation,”
Applied Numer-ical Mathematics , vol. 152, pp. 185–198, 2020.[12] M. I. Ahmad and P. Benner, “Interpolatory model reduction techniques forlinear second-order descriptor systems,” in
Proc. European Control Conf.ECC 2014, Strasbourg . IEEE, 2014, pp. 1075–1079.[13] S. Wyatt, “Issues in interpolatory model reduction: Inexact solves, secondorder systems and daes,” Ph.D. dissertation, Virginia Polytechnic Instituteand State University, Blacksburg, Virginia, USA, May 2012.1414] M. M. Rahman, M. M. Uddin, L. S. Andallah, and M. Uddin, “Tangentialinterpolatory projections for a class of second-order index-1 descriptor sys-tems and application to mechatronics,”
Production Engineering , pp. 1–11,2020.[15] P. Benner, J. Saak, and M. M. Uddin, “Balancing based model reductionfor structured index-2 unstable descriptor systems with application to flowcontrol,”
Numerical Algebra, Control and Optimization , vol. 6, no. 1, pp.1–20, 2016.[16] F. Tisseur and K. Meerbergen, “The quadratic eigenvalue problem,”
SIAMRev. , vol. 43, no. 2, pp. 235–286, 2001.[17] V. Mehrmann and T. Stykel, “Balanced truncation model reduction forlarge-scale systems in descriptor form,” 2005, chapter 20 (pages 357–361)of [19].[18] N. Truhar and K. Veseli´c, “Bounds on the trace of a solution to the Lya-punov equation with a general stable matrix,”
Syst. Cont. Lett. , vol. 56,no. 7–8, pp. 493–503, 2007.[19] P. Benner, V. Mehrmann, and D. C. Sorensen,