Frequency Limited \mathcal{H}_2 Optimal Model Reduction of Large-Scale Sparse Dynamical Systems
Xin Du, M. Monir Uddin, A. Mostakim Fony, Md. Tanzim Hossain, Mohammaed Sahadat-Hossain
aa r X i v : . [ m a t h . O C ] J a n Frequency Limited H Optimal Model Reductionof Large-Scale Sparse Dynamical Systems
Xin Du ∗ , M. Monir Uddin † , A. Mostakim Fony ‡ , Md. TanzimHossain § , and Mohammaed Sahadat-Hossain ¶ Abstract
We mainly consider the frequency limited H optimal model order re-duction of large-scale sparse generalized systems. For this purpose we needto solve two Sylvester equations. This paper proposes efficient algorithmto solve them efficiently. The ideas are also generalized to index-1 de-scriptor systems. Numerical experiments are carried out using PythonProgramming Language and the results are presented to demonstratethe approximation accuracy and computational efficiency of the proposedtechniques. keywords Frequency limited model reduction, H optimal condition,frequency limited Graminas, Sylvester and Lyapunov equations Model order reduction (MOR) is a process to approximate a high-order dy-namical system by a substantially low-order system with a maximum accuracy.This tool is now widely used in different disciplines of science, engineering andtechnology to reduce the complexity of the model. In general, the reducedorder models are used in controller design, simulation and optimization. Formotivation, applications and techniques of MOR see, e.g,. [1, 2].The commonly useful methods for the model reduction of large-scale lineartime invariant dynamical systems are the balanced truncation and H optimal ∗ 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 Mathematics and Physics, North south University, Dhaka-1229,Bangladesh, [email protected] H model reduction methods have been studied and investi-gated in [5, 6, 7, 8, 9]. See, also the references cited therein. In all these papersthe proposed technique is based on either Gramian assistance first-order opti-mality conditions [5] or the tangential interpolaton [6] of the transfer function.In fact both the conditions are coincided which is shown in [9]. These papersonly discuss the model reduction on the infinite frequency interval. For thetime limited case the we refer the readers to [10]. Although, in [11] authorsbriefly introduced optimal H model reduction problem of standard state spacesystems considering a restricted frequency interval, there the implementationdetails were not given.This paper focuses on the computational techniques of the frequency limitedoptimal H model reduction method of large-scale sparse systems. We mainlygeneralized the idea as in [12, 9] in which the proposed algorithm is called twosided iteration algorithm (TSIA). Moreover, to implement the frequency limitedTSIA we need to solve two frequency limited Sylvester equations. This paperalso discusses how to solve the Sylvester equations efficiently by preserving thesparsity of the system. Besides the generalized systems the idea is also extendedfor index-1 descriptor systems. The benefits of the algorithmic improvementspresented in this paper are illustrated by several numerical examples. We havegenerated the the results by using Python Programming Language. Rest of thispaper is organized as follows.Section 2 overview the TSIA and the optimal H model reduction of general-ized system. Then the ideas of this section are discussed in the next sections forthe frequency-limited model order reduction. Section 5 presents the algorithmfor solving frequency limited Sylvesters equations which provides projectors tocarry out the FLMOR. The results of the numerical experiments are depictedin Section 6 which show the efficiency and capability of the proposed methods. H optimal model order reduc-tion of generalized Systems The goal of this section is to review the basic idea of H optimal model or-der reduction of generalized Systems. Let us consider a linear time invariantcontinuous-time system of the form E ˙ x ( t ) = Ax ( t ) + Bu ( t ) ,y ( t ) = Cx ( t ) + Du ( t ) , (1)where E ∈ R n × n is non-singular, and A ∈ R n × n , B ∈ R n × p , C ∈ R m × n and D ∈ R m × p . The transfer-function matrix of this system is defined by G( s ) =2 ( sE − A ) − B + D , where s ∈ C . The controllability and the observabilityGramians of the system on the infinite frequency range can be defined as [13] P = 12 π Z ∞−∞ ( ıωE − A ) − BB T ( ıωE T − A T ) d ω and Q = 12 π Z ∞−∞ ( ıωE T − A T ) − C T C ( ıωE − A ) d ω, and they are the solutions of the continuous-time algebraic Lyapunov equations AP E T + EP A T + BB T = 0 and A T QE + E T QA + C T C = 0 , (2)respectively. We want to construct a substantially reduced-order model˙ˆ x ( t ) = ˆ A ˆ x ( t ) + ˆ Bu ( t ) , ˆ y ( t ) = ˆ C ˆ x ( t ) + ˆ Du ( t ) , (3)where ˆ A ∈ R r × r , ˆ B ∈ R r × p , ˆ C ∈ R m × r and ˆ D ∈ R m × p . The goal is to minimizethe error Ξ = k G − ˆG k H , (4)where k . k H denotes the system’s H -norm [1] and ˆ G ( s ) = ˆ C ( sI − ˆ A ) − ˆ B + ˆ D isthe transfer-function matrix of the reduced system. The H -norm of the errorsystem as defined in (4) can be measured byΞ =: q Tr (cid:0) C Ξ P Ξ C T Ξ (cid:1) or q Tr (cid:0) B T Ξ Q Ξ B Ξ (cid:1) , (5)where P Ξ and Q Ξ are the solutions of the Lyapunov equations A Ξ P Ξ E T Ξ + E Ξ P Ξ A T Ξ + B Ξ B T Ξ = 0 and (6a) A T Ξ Q Ξ E Ξ + E T Ξ Q Ξ A Ξ + C T Ξ C Ξ = 0 , (6b)in which E Ξ = (cid:20) E ˆ I (cid:21) , A Ξ = (cid:20) A ˆ A (cid:21) , B Ξ = (cid:20) B ˆ B (cid:21) and C Ξ = (cid:2) C − ˆ C (cid:3) . Now, partitioning P Ξ and Q Ξ as P Ξ = (cid:20) P MM T ˆ P (cid:21) and Q Ξ = (cid:20) Q NN T ˆ Q (cid:21) , respectivelyand plugging into (6) we obtain the following algebraic matrix equations AP E T + EP A T + BB T = 0 , (7a)ˆ A ˆ P + ˆ P ˆ A T + ˆ B ˆ B T = 0 , (7b) AM + EM ˆ A T + B ˆ B T = 0 , (7c) A T QE + E T QA + C T C = 0 , (7d)ˆ A T ˆ Q + ˆ Q ˆ A + ˆ C T ˆ C = 0 , (7e) A T N + E T N ˆ A − C T ˆ C = 0 , (7f)3here ˆ P and ˆ Q are respectively known as the controllability and observabilityGramians of the reduced systems. Therefore, the H norm of the error systemin (5) can be measured byΞ =: r Tr (
CP C T ) + Tr (cid:16) ˆ C ˆ P ˆ C T (cid:17) + 2Tr (cid:16) CM ˆ C T (cid:17) or r Tr ( B T QB ) + Tr (cid:16) ˆ B T ˆ Q ˆ B (cid:17) + 2Tr (cid:16) B T N ˆ B (cid:17) . (8)The first-order optimality conditions for the optimal H model reduction wasgiven in [14], which is known as Wilson conditions. Wilson conditions are basedon the first derivatives of (4) with respect to ˆ A , ˆ B and ˆ C as follows: ∇ Ξ ˆ A = 2( ˆ Q ˆ P + W T EV ) , ∇ Ξ ˆ B = 2( ˆ Q ˆ B + W T E − B ) , ∇ Ξ ˆ C = 2( ˆ C ˆ P − CV ) . Setting these three derivatives to zero leads to the
Wilson conditions ,ˆ Q ˆ P + N T EM = 0 , (9)ˆ Q ˆ B + N T E − B = 0 , (10)ˆ C ˆ P − CM = 0 . (11)These three conditions in fact yield the left and and right projection matrices tocompute an optimal reduced order system (3) and in the optimal reduced ordersystem the reduced matrices are formed asˆ A = W T E − AV, ˆ B = W T B, and ˆ C = CV, (12)where V = M ˆ P − and W T = − ˆ Q − N T and hence it can be proved that W T EV = I . However, we can not guarantee that ˆ P and ˆ Q are invertible,since to assure this the reduced model should be completely controllable andobservable [1]. In the case that they are invertible, the multiplication from theright is only a transformation of bases and does not change the subspace. Theidea by Xu and Zeng [9] was to satisfy the Wilson conditions by setting W = N and V = M. (13)Note that V and W can be computed by solving the matrix equations (7c) and(7f), respectively which are known as Sylvesters equations. Another importantobservation is that if we want to compute the optimal projection subspace wealready need the optimal solution ˆ A , ˆ B and ˆ C . However, this is not knownprior. A possible solution is to start with a reduced model, which emerged froman arbitrary projection of the original model, solve matrix equations (7c) and(7c), compute the projectors, and restart the process with the newly obtainedreduced model until we are satisfied. In this way we get a kind of a fixed pointiteration. This procedure is called two sided iteration algorithm (TSIA) by Xuand Zeng in [9].The Wilson conditions are Gramian-based conditions since they are relatedto Gramians of the system. The Hyland-Bernstein conditions [7] are anothergramian-based first-order optimal conditions, which were shown to be equiva-lent to the Wilson conditions [15]. Van Dooren et al., were characterized the4 lgorithm 1: Two-sided iteration algorithm (TSIA).
Input :
E, A, B, C, D . Output: ˆ A, ˆ B, ˆ C , ˆ D := D . Choose matrices W ∈ R n × r and V ∈ R n × r such that W T V = I . Construct the reduced-order matrices:ˆ A = W T E − AV , ˆ B = W T E − B and ˆ C = CV . while i ≤ N − do Compute V i and W i by solving Sylvesters AV + EV ˆ A T + B ˆ B T = 0 (14a) A T W + E T W ˆ A − C T ˆ C = 0 , (14b) Compute W i +1 = W i ( V Ti W i ) − and V i +1 = V i Construct the reduced-order matricesˆ A = W Ti +1 E − AV i +1 , ˆ B = W Ti +1 E − B and ˆ C = CV i +1 . i = i + 1. end while tangential interpolation based H optimal conditions in [12]. One drawback ofinterpolation based model reduction is to selection of interpolation points. How-ever in [15] authors proposed the Iterative Rational Krylov Algorithms (IRKA)to resolve this problem. On the other hand Xu and Zeng showed in [9] that boththe Gramian and interpolation based optimality conditions are the same. In [9]authors also presented the two sided iterative algorithm (TSIA) for a standardsystem which is slightly modified as in [12]. For our convenient the TSIA forthe H optimal model reduction for the generalized system (1) is summarizedin Algorithm 1 H optimal MOR of generalized systems This section moves to the frequency limited H optimal model reduction ofthe system (1). For this purpose we first define frequency limited Gramians.In the previous section the system Gramians have been defined on the infinitefrequency interval. If we replace the infinite interval (- ∞ , ∞ ) into a finiteinterval ω = [ ω , ω ]) then the controlability and the observablity Gramianscan be defined as P ω = 12 π Z ω ω ( ıνE − A ) − BB T ( ıνE T − A T ) d ν (15) Q ω = 12 π Z ω ω ( ıνE T − A T ) − C T C ( ıνE − A ) d ν, (16)which satisfies the frequency limited Lyapunov equations5 P ω E T + EP ω A T + B ω BB T + BB T B ∗ ω = 0 ,A T Q ω E + E T Q ω A + C ∗ ω C T C + C T CC ω = 0 , (17)with B ω = ı π ln( A + ıω E )( A + ıω E ) − and C ω = ı π ln( A + ıω E ) − ( A + ıω E ) . The goal in this paper is to constructed a reduced model ˆ G = ( ˆ A, ˆ B, ˆ C ) from agiven model G = ( E, A, B, C ) that minimizes the errorΞ ω = k G − ˆG k H ,ω , (18)where k . k H ,ω denotes the H -norm on the prescribed frequency range ω . Thetransfer-function matrix of the reduced system is The H -norm of the errorsystem as defined in (4) can be measured efficiently byΞ ω =: q Tr (cid:0) C Ξ P Ξ ,ω C T Ξ (cid:1) or q Tr (cid:0) B T Ξ Q Ξ ,ω B Ξ (cid:1) , (19)where P Ξ ,ω and Q Ξ ,ω are the solutions of the Lyapunov equations A Ξ P Ξ ,ω E T Ξ + E Ξ P Ξ ,ω A T Ξ + B Ξ ,ω B Ξ B T Ξ + B Ξ B T Ξ B ∗ Ξ ,ω = 0 and (20a) A T Ξ Q Ξ ,ω E Ξ + E T Ξ Q Ξ ,ω A Ξ + C ∗ Ξ ,ω C T Ξ C Ξ + C T Ξ C Ξ C Ξ ,ω = 0 , (20b)where B Ξ ,ω = ı π ln( A Ξ + ıω E Ξ )( A Ξ + ıω E Ξ ) − and C Ξ ,ω = ı π ln( A Ξ + ıω E Ξ ) − ( A Ξ + ıω E Ξ ) . Due to the structure of E Ξ , A Ξ , B Ξ and C Ξ we can Partition P Ξ ,ω , Q Ξ ,ω , B Ξ ,ω and C Ξ ,ω as follows P Ξ ,ω = (cid:20) P ω M ω M Tω ˆ P ω (cid:21) , Q Ξ ,ω = (cid:20) Q ω N ω N Tω ˆ Q ω (cid:21) ,B Ξ ,ω = (cid:20) B ω
00 ˆ B ω (cid:21) , C Ξ ,ω = (cid:20) C ω
00 ˆ C ω (cid:21) . Therefore (20) yields AP ω E T + EP ω A T + B ω BB T + BB T B ∗ ω = 0 , (21a)ˆ A ˆ P ω + ˆ P ω ˆ A T + ˆ B ω ˆ B ˆ B T + ˆ B ˆ B T ˆ B ∗ ω = 0 , (21b) AM ω + EM ω ˆ A T + B ω B ˆ B T + B ˆ B T ˆ B ∗ ω = 0 , (21c) AQ ω E T + EQ ω A T + C ∗ ω C T C + C T CC ω = 0 , (21d)ˆ A ˆ Q ω + ˆ Q ω ˆ A T + ˆ C ∗ ω ˆ C T ˆ C + ˆ C T ˆ C ˆ C ω = 0 , (21e) AN ω + EN ω ˆ A T + C ∗ ω C T ˆ C − C T ˆ C ˆ C ω = 0 , (21f)6 lgorithm 2: Two-sided iteration algorithm (TSIA).
Input :
E, A, B, C . Output: ˆ A, ˆ B, ˆ C , ˆ D := D . Choose matrices W ∈ R n × r and V ∈ R n × r such that W T V = I . Construct the reduced-order matrices ˆ A = W T E − AV , ˆ B = W T E − B and ˆ C = CV . while ( i ≤ N − ) do Compute V i = M ω and W i = N ω by solving the Sylvester AM ω + EM ω ˆ A T + B ω B ˆ B T + B ˆ B T ˆ B ∗ ω = 0 , (22a) AN ω + EN ω ˆ A T + C ∗ ω C T ˆ C − C T ˆ C ˆ C ω = 0 , (22b)Compute W i +1 = W i ( V Ti W i ) − and V i +1 = V i . Construct the reduced-order matricesˆ A = W Ti +1 E − AV i +1 , ˆ B = W Ti +1 E − B and ˆ C = CV i +1 . i = i + 1. end while with ˆ B ω = ı π ln( ˆ A + ıω I )( ˆ A + ıω I ) − andˆ C ω = ı π ln( ˆ A + ıω I ) − ( ˆ A + ıω I ) . Following the discussion in the above section here we also construct reducedorder model by constructing the reduced matrices as in (12). We solve the theSylvesters equations (21c) and (21f) to construct V = M ω and W = N ω . Theconstructed reduced order model is H optimal on the limited frequency rangeand satisfies Wilson’s first-order optimality conditions. The whole procedure issummarized in Algorithm 2. The main computation tasks in this algorithm isto solve sparse-dense Sylvesters equations (22a) and (22b). Following sectionwill presents how to solve them efficiently. Above section shows that to perform the frequency limited model reduction ofsystem (1) we need to solve two frequency limited matrix equations namelySylvester equations (22a) and (22b). This section discusses how to solve themefficiently. Since the Sylvester equations (22a) and (22b) are duel of each otherwe only interested to elaborate the solution of (22a). Another one can be solvedby applying the same procedure. For our convenient we rewrite the equation(22a) as AX + EX ˆ A T + F = 0 , (23)where F = B ω BB T + BB T B ∗ ω and X = M ω . The technique that we havefollowed here was presented in [16] where E = I is an identity matrix. In [17]7 lgorithm 3: Solution of semi-generalized Sylvester equations.
Input :
E, A, ˆ A, F from (23).
Output: X ∈ R n × r solution of (23). Compute the Schur decomposition ˆ A = QSQ ∗ and Define ˜ F = F Q for i = 1 , · · · , r do Compute ˆ F = − ˜ F j − E P j − i =1 S i,j ˜ X i Solve ( A + S jj E ) ˜ X j = ˆ F end for X = ˜ XZ ∗ .authors generalized the idea of [16] for the equation like (23) where F = B ˆ B .Considering the Schur decomposition of ˆ A as QSQ ∗ such that QQ ∗ = Q ∗ Q = I and inserting this into (23) we get AX + EXQSQ ∗ + F = 0 . (24)By multiplying this from the right with Q we obtain A XQ |{z} ˜ X + E XQ |{z} ˜ X S + F Q |{z} ˜ F = 0 , (25)Observing that S is a upper triangular matrices leads to a formula for the firstcolumn of ˜ X : A ˜ X + E ˜ X S , + ˜ F = 0 , (26) ⇔ ( A + S E ) ˜ X = − ˜ F . (27)For all other columns we have to take care of the linear combination of E matrix.If we consider the second column of the solution( A + S E ) ˜ X = − ˜ F − S E ˜ X . (28)In this way the arbitrary column j of ˜ X we find( A + S jj E ) ˜ X j = − ˜ F j − E j − X i =1 S ij ˜ X i . (29)To obtain the solution of the original system we multiply ˜ X by U ∗ from theright. The index 1 descriptor system that we consider in the section has the followingform. E ˙ x ( t ) = J x ( t ) + J z ( t ) + B u ( t )0 = J x ( t ) + J z ( t ) + B u ( t ) y ( t ) = C x ( t ) + C z ( t ) + D a u ( t ) , (30)8here x ( t ) ∈ R n is the vector with differential variables and z ( t ) ∈ R n isthe vector with algebraic variables. Model reduction of such descriptor systemhas been discussed in a couple of previous research articles, e.g., [18, 19, 20, 21,22, 20, 23, 24] on an unrestricted frequency limit. In [18] author uses spectralprojector to split the system into finite and infinite parts and balancing basedmethod was applied to the finite part. Other papers implemented MOR with-out computing the spectral projector, rather eliminating the algebraic part thesystem was converted into an ODE system. However, practical implementa-tion was carried out without computing the ODE system explicitly. This papergeneralizes this idea for of the FLMOR discussed in the above section.By eliminating the algebraic variables i.e., z ( t ) ∈ R n of the system we obtaina generalized system (1) where the coefficients matrices are defined as E := E , A := J − J J − J , B := B − J J − B ,C := C − C J − J , D := D a − C J − B . (31)The index-1 and generalized systems are equivalent since the responses of thesystems are same and their finite eigenvalues are coincided. Such structured sys-tem are arising in power system model [25]. For the FLMOR of index-1 systemwe define V and W by solving the corresponding Sylvesters equations of the gen-eralized system as discussed in Section 3. Now, applying these transformationsthe reduced system matrices can be constructed as:ˆ A := ˆ J − ˆ J J − ˆ J , ˆ B := ˆ B − ˆ J ˆ J − B , ˆ C := ˆ C − C J − ˆ J , ˆ D := D a − C J − B , (32)where ˆ J = W T E − J V , ˆ J = W T J , ˆ J = J V , ˆ B = W T E − B andˆ C = C V . To compute the the projectors V and W by solving the corre-sponding Sylvesters equations is a challenging task since the input matrices in(32) are highly dense. In the following text we discuss how to solve the Sylvestersequations related to the index-1 system (30) efficiently.To solve the Sylvesters equations we can use Algorithm 3. In this algorithmthe main expensive task is to solve a linear system at each iteration step. AtStep 4 of the algorithm we need to solve the linear system( A + S ii E ) ˜ X j = ˆ F .
Plugging A and E from (31) we obtain( J − J J − J + S ii E ) ˜ X j = ˆ F , which can be rewritten as( J + S ii E − J J − J ) ˜ X j = ˆ F . (33)A close observation reveal that instead of this we can solve the following linearsystem (cid:20) J + S ii E J J J (cid:21) (cid:20) ˜ X j Γ (cid:21) = (cid:20) ˆ F (cid:21) (34)for ˜ X j . Although, linear system in (34) is larger than the system in (33) it issparse and hence can be solved by any sparse solvers (e.g., direct [26, 27] oriterative e.g., [28, 29]) efficiently. 9 lgorithm 4: Sylvester equation for index-1 system.
Input : E , J , J , J , J , B , B , ˆ A . Output: X ∈ R n × r . Form F = B ω B ˆ B T + B ˆ B T ˆ B ∗ ω , where B ω = ı π ln( J + ıω E − J J − J )( J + ıω E − J J − J ) − , ˆ B ω = ı π ln( ˆ A + ıω I ) − ( ˆ A + ıω I ) , and B = B − J J − B Compute the Schur decomposition ˆ A = QSQ ∗ and Define ˜ F = F Q for i = 1 , · · · , r do Compute ˆ F = − ˜ F j − E P j − i =1 S i,j ˜ X i Solve (cid:20) J + S ii E J J J (cid:21) (cid:20) ˜ X j Γ (cid:21) = (cid:20) ˆ F (cid:21) for ˜ X j end for X = ˜ XZ ∗ . To asses the efficiency of our proposed techniques in this section we discuss thenumerical results. For our convenient we have splitted the section into severalsubsections.
We consider the following model examples for the numerical experiments.
Example 6.1 (International Space Station (ISS)) . This is a model of stage 1R(Russian Service Module) of the ISS. It has n=270 states, p=3 inputs and m=3outputs. The details of the model can be obtain in [30].
Example 6.2 (Clamped beam model (CBM)) . This structural model was ob-tained by spatial discretization of an appropriate partial differential equation(see [31]). The dimension of the model is n=348 and it is a single input singleoutput i.e., SISO system. The input represents the force applied to the structureat the free end, and the output is the resulting displacement.
Example 6.3 (Triple chain oscillator (TCO) model) . Although this examplewas originated in [32], the setup was described in [33] which resulted in second-order system. We convert into firs-order form in which the dimension of thesystem is n=10000. It is also an SISO system and input, output matrices aretranspose of each other.
Example 6.4 (Power system model) . We consider several Brazilian Power Sys-tem (BPS) models from [21] which are in index-1 descriptor form. Table 1 shows10 odel name n n m/p BPS-606 606 1142BPS-1142 1142 8593 4/4BPS-1450 1450 9855BPS-1693 1693 b11582Table 1: Dimension of differential and algebraic variables of different Brazilianpower system models.Model r Ξ ω ΞISS 30 7.3244 × − × − TCO 30 3.4850 × − × − BPS-606 30 7.8992 × − × − BPS-1142 35 3.4421 × − × − BPS-1450 35 9.3818 × − × − BPS-1693 45 4.2995 × − × − Table 2: Dimension of reduced order models and H norms of the error systemswith and without limited frequency intervals.number of differential ( n ) and algebraic n variables and inputs/outputs of sev-eral models which are used for the numerical experiments here. The experiments are carried out with Python 3.7.9 on a board with AMDRyzen TM Threadripper TM The proposed techniques was applied to the all model examples mentionedabove. For the ISS, CBM, TCO we apply Algorithm 2 to obtain frequencylimited reduced order model. On the other for the BPS models we have appliedthe techniques discussed in Section 5. We have computed different dimensionalreduced order models for different model examples which are mentioned in Ta-ble 2. The table also shows H norm of the error systems in both frequencyrestricted and unrestricted cases. For all the models frequency restricted re-duced order model show much more better accuracy than the frequency fre-quency underused ones within the assigned frequency intervals. We also haveinvestigated the frequency domain analysis including the errors of the originaland reduced order models by using sigma plots. Exemplary, we have depictedthe frequency responses of TCO and BPS models only. Figures 1 and 2 showcomparisons of the frequency responses of the of the original and the reducedorder models of TCO and BPS models, respectively . In both the figures fromabsolute and the relative errors we observe that frequency restricted reducedorder models approximate with the original models with higher accuracy withinon the prescribed frequency ranges. 11 0 . . . . ω σ m a x ( G ( j ω )) Full model Frequency unrestricted Frequency restricted (a) Sigma plot . . . . − ω σ m a x ( G ( j ω ) - ˆ G ( j ω )) (b) Absolute error . . . . − − ω σ m a x ( G ( j ω ) - ˆ G ( j ω )) σ m a x ( G ( j ω )) (c) Relative error Figure 1: Comparison of original and 30 dimensional reduced systems on thefrequency range [1,2] for TCO. 12 4 6 8 10 12 14 1610 . ω σ m a x ( G ( j ω )) Full model Frequency unrestricted Frequency restricted (a) Sigma plot − − ω σ m a x ( G ( j ω ) - ˆ G ( j ω )) (b) Absolute error − − ω σ m a x ( G ( j ω ) - ˆ G ( j ω )) σ m a x ( G ( j ω )) (c) Relative error Figure 2: Comparison of original and 45 dimensional reduced systems on thefrequency range [6,10] for BPS-1693. 13PS-606 BPS-1142 BPS-1450 BPS-169302468 0.62 1.94 3.19 7.740.53 0.73 0.87 1.16 t i m e ( s ec ) Algorithm 3Algorithm 4Figure 3: Comparison of computational time to solve the Sylvester equation fordifferent dimensional BPS model using Algorithms 3 and 4.
We already mentioned that in the TSIA the main computational task is solvingtwo Sylvester equations. This paper presents Algorithm 4 to solve the Sylvesterequations efficiently for index-1 system. We know that by converting index-1 system into a generalized system we can solve the Sylvester equation usingAlgorithm 3. Figure 3 shows the time comparisons of Algorithms 3 and 4 forsolving Sylvester equations of index-1 system. We see that if the dimensionof the system is increased then the computational times for Algorithms 3 isincreased rapidly. On the other hand the computational time with Algorithm 4is nominal in compare to Algorithms 3.
In this paper we have discussed the frequency limited H optimal model orderreduction of large-scale sparse dynamical systems. For this purpose we havesolved two mixed (Sparse and dense) Sylvester equations which are dual of eachother. We have shown that how to solve the Sylvester equations efficientlywithout loosing the sparsity of large sparse system matrices. The ideas are alsogeneralized to index-1 descriptor systems. Index-1 system can be converted intoa generalized system by eliminating the algebraic equations, which however con-vert the system from sparse to dense. We have discussed how to perform themodel reduction without converting the dense system explicitly. Numerical ex-periments has been carried out to demonstrate the approximation accuracy andcomputational efficiency of the proposed algorithm using Python ProgrammingLanguage. 14 cknowledgments. This research work was funded by NSU-CTRG researchgrant under the project No.: CTRG-19/SEPS/05. It was also supported byNational Natural Science Foundation of China under Grant No. (61873336,61873335), the Fundamental Research Funds for the Central Universities underGrant (FRF-BD-19-002A), and the High-end foreign expert program of Shang-hai University,
References [1] A. Antoulas,
Approximation of Large-Scale Dynamical Systems , ser. Ad-vances in Design and Control. Philadelphia, PA: SIAM Publications,2005, vol. 6.[2] M. M. Uddin,
Computational Methods for Approximation of Large-ScaleDynamical Systems . New York, USA: Chapman and Hall/CRC, 2019.[3] W. Gawronski and J. Juang, “Model reduction in limited time and fre-quency intervals,”
Int. J. Syst. Sci. , vol. 21, no. 2, pp. 349–376, 1990.[4] P. Benner, P. K¨urschner, and J. Saak, “Frequency-limited balanced trunca-tion with low-rank approximations,”
SIAM J. Sci. Comput. , vol. 38, no. 1,pp. A471–A499, Feb. 2016.[5] P. Van Dooren, K. Gallivan, and P.-A. Absil, “ H -optimal model reductionof MIMO systems,” Appl. Math. Lett. , vol. 21, pp. 1267–1273, 2008.[6] 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.[7] D. Hyland and D. Bernstein, “The optimal projection equations for modelreduction and the relationships among the methods of wilson, skelton, andmoore.”[8] W. Y. Yan and J. Lam, “An approximate approach to H optimal modelreduction,” IEEE Trans. Autom. Control , vol. 44, no. 7, pp. 1341–1358,1999.[9] Y. Xu and T. Zeng, “Optimal H model reduction for large scale MIMOsystems via tangential interpolation,” International Journal of NumericalAnalysis and Modeling , vol. 8, no. 1, pp. 174–188, 2011.[10] P. Goyal and M. Redmann, “Time-limited H -optimal model order reduc-tion,” Applied Mathematics and Computation , vol. 355, pp. 184–197, 2019.[11] D. Petersson and J. L¨ofberg, “Model reduction using a frequency-limited H -cost,” Systems & Control Letters , vol. 67, pp. 32–39, 2014.[12] P. Van Dooren, K. A. Gallivan, and P. A. Absil, “H2-optimal model reduc-tion of mimo systems,”
Appl. Math. Lett. , vol. 21, no. 12, pp. 1267–1273,2008. 1513] K. Zhou, G. Salomon, and E. Wu, “Balanced realization and model reduc-tion for unstable systems,”
Internat. J. Robust and Nonlinear Cont. , vol. 9,no. 3, pp. 183–198, 1999.[14] D. A. Wilson, “Optimum solution of model-reduction problem,”
Proceed-ings of the Institution of Electrical Engineers , vol. 117, no. 6, pp. 1161–1165,1970.[15] S. Gugercin, A. C. Antoulas, and C. Beattie, “H 2 model reduction forlarge-scale linear dynamical systems,”
SIAM journal on matrix analysisand applications , vol. 30, no. 2, pp. 609–638, 2008.[16] D. C. Sorensen and A. C. Antoulas, “The sylvester equation and approxi-mate balanced reduction,”
Numer. Lin. Alg. Appl. , vol. 351–352, pp. 671–700, 2002.[17] P. Benner, M. K¨ohler, and J. Saak, “Sparse-dense sylvester equations in h -model order reduction,” 2011.[18] T. Stykel, “Analysis and numerical solution of generalized Lyapunov equa-tions,” Ph.D. Thesis, Technische Universit¨at Berlin, Berlin, 2002.[19] 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.[20] M. M. Uddin, J. Saak, B. Kranz, and P. Benner, “Computation of a com-pact state space model for an adaptive spindle head configuration with piezoactuators using balanced truncation,”
Production Engineering , vol. 6, pp.577–586, 2012.[21] F. Freitas, J. Rommes, and N. Martins, “Gramian-based reduction methodapplied to large sparse power system descriptor models,”
IEEE Trans.Power Syst.
Numerical Algebra, Control & Optimiza-tion , vol. 9, no. 2, p. 173, 2019.[24] M. M. Uddin, M. S. Hossain, and M. F. Uddin, “Rational krylov subspacemethod (rksm) for solving the lyapunov equations of index-1 descriptorsystems and application to balancing based model reduction,” in .IEEE, 2016, pp. 451–454.[25] R. W. Freund, “Structure-preserving model order reduction of RCL circuitequations,” in
Model Order Reduction: Theory, Research Aspects and Ap-plications , W. H. A. Schilders, H. A. van der Vorst, and J. Rommes, Eds.Springer-Verlag, 2008, pp. 49–73. 1626] T. A. Davis,
Direct Methods for Sparse Linear Systems , ser. Fundamentalsof Algorithms. Philadelphia, PA, USA: SIAM, 2006, no. 2.[27] I. S. Duff, A. M. Erisman, and J. K. Reid,
Direct methods for sparse ma-trices . Oxford, UK: Clarendon Press, 1989.[28] H. A. Van der Vorst,
Iterative Krylov Methods for Large Linear Systems .Cambridge: Cambridge University Press, 2003.[29] Y. Saad,
Iterative Methods for Sparse Linear Systems . Philadelphia, PA,USA: SIAM, 2003.[30] S. Gugercin, A. C. Antoulas, and N. Bedrossian, “Approximation of theinternational space station 1r and 12a models,” in
Proc. of the 40th IEEEConferences on Decision and Control , 2001, pp. 1515–1516.[31] A. C. Antoulas, D. C. Sorensen, and S. Gugercin, “A survey of modelreduction methods for large-scale systems,”
Contemp. Math. , vol. 280, pp.193–219, 2001.[32] N. Truhar and K. Veseli´c, “An efficient method for estimating the optimaldampers’ viscosity for linear vibrating systems using Lyapunov equation,”