Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Parametric model order reduction of thermal models using the bilinear interpolatory rational Krylov algorithm

Parametric model order reduction of thermal models using the bilinear interpolatory rational... Mathematical and Computer Modelling of Dynamical Systems, 2015 Vol. 21, No. 2, 103–129, http://dx.doi.org/10.1080/13873954.2014.924534 Parametric model order reduction of thermal models using the bilinear interpolatory rational Krylov algorithm a b Angelika Bruns * and Peter Benner a b Robert Bosch GmbH, CR/ARH2, Robert-Bosch-Platz 1, 70839 Gerlingen, Germany; Max Planck Institute for Dynamics of Complex Technical Systems, Computational Methods in Systems and Control Theory, Sandtorstr. 1, 39106 Magdeburg, Germany (Received 25 September 2013; accepted 12 May 2014) The Bilinear Interpolatory Rational Krylov Algorithm (BIRKA; P. Benner and T. Breiten, Interpolation-based H -model reduction of bilinear control systems, SIAM J. Matrix Anal. Appl. 33 (2012), pp. 859–885. doi:10.1137/110836742) is a recently developed method for Model Order Reduction (MOR) of bilinear systems. Here, it is used and further developed for a certain class of parametric systems. As BIRKA does not preserve stability, two different approaches generating stable reduced models are pre- sented. In addition, the convergence for a modified version of BIRKA for large systems is analysed and a method for detecting divergence possibly resulting from this modifica- tion is proposed. The behaviour of the algorithm is analysed using a finite element model for the thermal analysis of an electrical motor. The reduction of two different motor models, incorporating seven and thirteen different physical parameters, is performed. Keywords: parametric model order reduction; stability preservation; thermal heat transfer; finite element modelling Subject Classification: 34K17; 93A15 1. Introduction The design of a new product in industry is often accompanied by computer simulations, which are an excellent tool for analysing the behaviour of the involved components. For the estimation of the product’s lifespan, thermal simulations are employed to detect whether an excessive increase in temperature will lead to a fatigue or even deterioration of the materials. Thermal simulations are often carried out using the finite element method, for which excellent tools exist. After a finite element discretization, the tool solves a system of differential equations of the following form: Ext _ðÞ ¼ AxðÞ t þ BuðÞ t ; (1) ytðÞ ¼ CxðÞ t ; n nn nm pn with the state vector x(t) 2 R , the constant matrices E; A 2 R , B 2 R , C 2 R , m p the inputs uðtÞ2 R , and outputs yðtÞ2 R . Usually, the thermal finite element models are complex (often more than 500,000 degrees of freedom), and require therefore a large numerical effort to converge to a solution. *Corresponding author. Email: angelika.bruns@de.bosch.com © 2014 Taylor & Francis 104 A. Bruns and P. Benner Model Order Reduction (MOR) is a powerful method to reduce the dimension of the underlying system and therefore the simulation time significantly. First, an approximation nr of the original state vector xtðÞ  Vx ðÞ t with V 2 R is required while the dimension r of x (t) should be much smaller than the original dimension n. Equation (1) now reduces to EVx_ ðÞ t ¼ AVx ðÞ t þ BuðÞ t þ "ðÞ t ; r r (2) y ðÞ t ¼ CVx ðÞ t ; r r with a residual ε(t). By multiplication from the left with the transpose of a matrix nr W 2 R which is chosen such that the residual ε(t) is orthogonal to the subspace spanned by its columns (i.e. W ε(t) = 0), the following reduced system is obtained: E A B r r r zfflfflffl}|fflfflffl{ zfflfflffl}|fflfflffl{ zffl}|ffl{ T T T W EV x_ ðÞ t ¼ W AV x ðÞ t þ W B utðÞ; r r (3) y ðÞ t ¼ CV x ðÞ t ; r r |{z} rr rm pr p1 where E , A 2 R , B 2 R , C 2 R and y ðtÞ2 R . The goal of projection- r r r r r based MOR is to determine the matrices V and W. For this class of systems, there exist many well-known methods such as Balanced Truncation, Proper Orthogonal Decomposition (POD) or Krylov subspace methods. An overview of these methods can be found in [1]. In the context of thermal simulations, several physical effects have to be considered. Heat sources, convection between the component and the surrounding environment and predefined temperatures on certain parts have to be taken into account. For this reason, the system (1) has to be reformulated to include the parameters h and the loads u(t) corresponding to the considered physical effects above: Ext _ðÞ ¼ A þ h N xtðÞþ BuðÞ t ; i i (4) i¼1 ytðÞ ¼ CxðÞ t ; nn nl kn n l k with E, A, N 2 R , B 2 R , C 2 R , xðtÞ2 R , uðtÞ2 R and yðtÞ2 R . Different parameters and loads denote different situations and designs. It is therefore desirable to reduce the system while keeping the dependence on the parameters h . This problem belongs to the class of parametric systems, for which reduction – the so-called parametric Model Order Reduction (pMOR) – was examined by a number of authors. Using a Moment-Matching approach, Weile et al. [2] reduced systems with two parameters, while Daniel et al. [3] extended this method to the multiple parameter case. Another well-known idea is the calculation of different reduced models for different parameter values followed by an interpolation between their underlying matrices. For example, Baur et al. [4] combine this idea with interpolatory MOR, Panzer et al. [5] use Krylov subspace methods and Amsallem et al. develop an interpolation procedure on a certain mani- fold [6]. Recently, Breiten and the second author proposed to rewrite a parametric system (4) as a bilinear system [7] of the form: Mathematical and Computer Modelling of Dynamical Systems 105 Ext _ðÞ ¼ AxðÞ t þ N u xtðÞþ BuðÞ t ; i i i¼1 (5) bil ytðÞ ¼ CxðÞ t ; Instead of reducing a linear, parametric system, the reduction of a bilinear system with parameters as inputs can now be realized. Al-Baiyat and Bettayeb studied the reduction of bilinear systems based on balancing, see for example [8], whereas Bai and Skoogh developed methods for reducing them with Krylov subspace methods [9]. This approach was further examined by Breiten and Damm [10]. An H optimal method was developed by Zhang and Lam [11] and transferred to interpolatory MOR by Breiten and the second author [12]. This new algorithm, the so-called Bilinear Interpolatory Rational Krylov Algorithm (BIRKA), will be examined within an industrial context in this work and will be applied to the reduction of an electrical motor. Other methods for the efficient simulation of electrical motor models can be found in the following works: Bertin et al. [13] and Gao et al. [14] apply their reduction methods to thermal networks. As it is done with BIRKA in this work, thermal finite element models of a motor have been reduced by Zhou et al. [15] by capturing the most significant eigenmodes and only considering the temperature ‘hot spots’. The outline of the paper is as follows: In Section 2, BIRKA is introduced and its operation mode for large systems is explained. A norm estimate which can detect divergence in certain cases is introduced in Section 3, while two approaches for preser- ving the system’s stability are introduced in Section 4, as stability is not preserved by the algorithm. Section 5 introduces the model – an electrical motor, and the underlying thermal modelling, while Section 6 shows how the model is parameterized. Results of the reduction with BIRKA can be found in Section 7. 1.1. Notation and definitions nm kl The Kronecker product of two matrices A 2 R and B 2 R is defined as 0 1 a B ... a B 11 1m B . . C . . A#B ¼ : @ A . . aBa B n1 nm nm The vectorization operator of a matrix A 2 R is: vecðÞ A ¼½ a ... a a ... a ... a ... a 11 n1 12 n2 1m nm nn The notation λ (A),i= 1, …, n, is used for the eigenvalues of the matrix A 2 R and λ i i nn (A, B),i= 1, …, n, for the eigenvalues of the pencil (A, B) with A; B 2 R . Similarly, nm nn σ (A),i= 1, …, n, are the singular values of the matrix A 2 R . If the matrix A 2 R is positive (semi)definite, this is written as A ≥ 0or A> 0, respectively. The analogue notation holds for negative (semi)definite matrices. The spectral norm of a matrix is denoted as kk M ¼ σ ðÞ M , and the spectral norm condition number of a matrix as max σ ðÞ M 1 max κ ðÞ M ¼kk M M ¼ , where σ ðÞ M =σ ðÞ M denotes the maximal/ 2 max min 2 σ ðÞ M 2 min minimal singular value of M. 106 A. Bruns and P. Benner Definition 1.1. For a linear system (1) the H -norm is defined as: 2 T kk H ¼ tr HðÞ iω HiðÞ ω dω (6) 2π −1 where H(s)= C(sE − A) B is the transfer function corresponding to the linear system (1). Definition 1.2. For a bilinear system (5) the H -norm is defined as: Z Z 1 m 1 1 X X ðÞ l ;...l ðÞ l ;...l 2 1 k 1 k ¼ tr ... g ðg Þ ds ... ds (7) kk bil 1 k H k k 0 0 k¼1 l ;...; l ¼1 1 k As As ðÞ l ;...l 1 1 k As N e k1 N ...e k l l 1 2 with g ðÞ s ; ... ; s¼ Ce b : 1 k l k k 2. The bilinear interpolatory rational Krylov algorithm The objective of MOR is to generate a reduced model that approximates the original one with as small deviation as possible. This deviation can be measured, for example, in terms of the error in the H -norm. For a bilinear system (5), necessary H optimality conditions, 2 2 which minimize this error locally, have been derived in [11,12]. The reduced models generated with BIRKA fulfil these conditions. The reader should note that for linear, parametric systems optimality conditions are not completely known yet. A first result can be found in [4], where the condition is given using a norm on H #L ðÞ Ω , with parameter 2 2 domain Ω. The finite element discretization of industrial models leads to a system with E ≠ I . Due to their large dimension, the inversion of E would be numerically −1 expensive or even impossible and would in addition lead to full matrices E A and −1 E B. Hence, models with nonsingular E ≠ I are considered, but explicit inversion of E needs to be avoided. This can be accomplished by a modification of the original BIRKA [12] which incorporates the E in the iteration as it can be seen in Algorithm 1. There, projection matrices V and W are iteratively computed. A set of starting matrices (E = I,A ,(N ) ,B,C ) of the reduced order r is required, which is chosen randomly or r r r r k r r by an informed guess. The stopping criterion is based on the change of the eigenvalues Λ −1 of the diagonalization A = S ΛS of the iteratively generated matrices. The algorithm is stopped if the change in the eigenvalues Λ of two subsequent models is smaller than a predefined tolerance. The projection matrices for the generation of the next reduced model are calculated by T T e ~ vecðÞ V ¼ I # A  Λ# E  N # N B # B vecðÞ I ; (8a) r k k m k¼1 Mathematical and Computer Modelling of Dynamical Systems 107 T T T T e ~ vecðÞ W ¼ I # A  Λ# E  N # N C # C vec I : (8b) r k k p k¼1 Algorithm 1 [12] Bilinear IRKA for systems with E ≠ I, E nonsingular Input: E, A, N ,B,C,E,A,(N ) ,B,C k r r r k r r opt opt opt opt opt Output: E ¼ I ; A ;ðÞ N ; B ; C r r r r k r r 1: while (Change in Λ >0) do 1 T 1 1 ~ ~ ~ 2: A ¼ SΛS ; B ¼ S B ; C ¼ C SN ¼ S ðÞ N S r r r k r T T ~ ~ 3: vecðÞ V ¼ I #A  Λ#E  N #N B #B vecðÞ I r k k m k¼1 1 T T T T T 4: vecðÞ W ¼ I #A  Λ#E  N #N C #C vec I r k k p k¼1 5: V = orth(V), W = orth(W) % orth computes an orthonormal basis 1 1 1 T T T T T T 6: A ¼ðÞ W EV W AV ; ðÞ N ¼ðÞ W EV W N V ; B ¼ðÞ W EV W B; C ¼ CV r r k r r 7: end while opt opt opt opt 8: A ¼ A ;ðÞ N ¼ðÞ N ; B ¼ B ; C ¼ C r r r r r r k k r r Due to the Kronecker product, the calculation of the projection matrices V and W is not possible for large systems. In [12] an iterative method was proposed to overcome this difficulty. For the derivation a Neumann Series is used in the following way: T T ~ ~ vecðÞ V ¼ I # A  Λ# E  N # N B # B vecðÞ I r k k m k¼1 "# ! 1 m X X ¼ ðÞ I # A  Λ# E N # N r k k |{z} i¼0 k¼1 ðÞ ðÞ I # A  Λ# E B # B vecðÞ I r m ¼ðÞ I # A  Λ# E B # B vecðÞ I þ r m |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V (9) ðÞ I # A  Λ# E N # N vec V r k k k¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V j1 ...þðÞ I # A  Λ# E N # N vec V þ ... r k k k¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V ¼ vec V ; j¼1 108 A. Bruns and P. Benner 1 m where (×) is only valid if ðÞ I # A  Λ# E ð N # N Þ <1 holds. In prac- r k k k¼1 tice, the infinite sum is truncated after an appropriate number of additions. The columns of the summands V are now calculated without using any Kronecker products: V ¼ðÞ λ E  A BB i i 2 1 V ¼ðÞ λ E  A N V N i k k k¼1 j 1 j1 V ¼ðÞ λ E  A N V N ; for i ¼ 1; ... ; r: i k k k¼1 This calculation can be executed in the same way for vec(W). The same projection matrices are calculated using the Truncated BIRKA proposed by Flagg [16]. The large matrices (−λ E − A) can be factorized by an LU-decomposition so that V can be efficiently calculated. For large models (n >10 ), this approach might not be feasible due to its memory demands. One can then solve the linear systems by an iterative method (cf. Saad [17]). 3. Estimation of the norm for the Kronecker product approximation Approximating the Kronecker product as in (9) can lead to divergence if ðÞ I # A  Λ# E N # N   1: r k k k¼1 It is advisable to check whether this norm remains smaller than one during the execution of BIRKA, as divergence might lead to poor reduced order models. However, a direct rnrn calculation of the norm would involve the inversion of ðI #A  Λ#EÞ2 R which might not be feasible for large systems due to memory problems. Hence, the calculation of the Kronecker product has to be avoided. To this aim, the following norm estimation is introduced: ðÞ I # A  Λ# E ð N # N Þ r k k k¼1 ðÞ I # A  Λ# E ð N # N Þ r k k k¼1 (10) ðÞ I # A  Λ# E  N # N r k k k¼1 ðÞ I # A  Λ# E N kk N : r k k |{z} k¼1 see below If the last expression is smaller than 1, the algorithm is definitely usable. We have thus derived a sufficient condition. Mathematical and Computer Modelling of Dynamical Systems 109 The normðÞ I # A  Λ# E can be calculated without explicit inversion of the matrix. The following Lemmata (cf [18]. Chapter I.4 and [1] Chapter 3) will be used to establish the new result for the calculation of the corresponding norm in Proposition 3.3. nn Lemma 3.1. For M 2 C nonsingular, M ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : min λ MM i¼1...n i Lemma 3.2. For a normal matrix M: M ¼ : min jj λðÞ M i¼1...n i By using the two Lemmata 3.1 and 3.2, we derive the following proposition, which will be used for the calculation of the norm ofðÞ I # A  Λ# E . nn Proposition 3.3. For A; E 2 R , symmetric, D = diag(d , …,d ): 1 r ðÞ I # A  D# E ¼ ; where jλ ðA  d EÞj min k for d 2 R qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi θ ¼ min : for d 2 C k¼1...r λ ððA  d EÞðA  d EÞ Þ k min k k Proof. The above matrix can be written as follows: 2 3 A  d E 0 ... 0 6 7 0 A  d E ... 0 6 7 ðÞ I #A  D#E¼ 6 . . 7; . . 4 5 . . 0 ... 0 A  d E with d 2 C.For d 2 R it is obvious that (−A − d E) is normal due to A, E symmetric, k k and thus Remark 3.2 can be used for calculating λ (−A − d E). For d 2 C the min k eigenvalue λ ðÞ A  d EðÞ A  d E is determined using Lemma 3.1. Taking the min k k minimum of all calculated eigenvalues and inverting it concludes the proof. ⊚ The calculation ofðÞ I # A  Λ# E can now be done by Proposition 3.3 using the MATLAB [19] function eigs. For the estimation of the norm as given in (10) it remains to calculate the norms of N and N , which is done in MATLAB with the functions k k normest and norm, respectively. For randomly chosen initial values the norm estimate is possibly greater than 1. However, as Λ and N change towards their optimal values, the norm estimate improves. At least two or three iterations should therefore be performed to check whether the norm for better approximations of Λ and N is smaller than 1. k 110 A. Bruns and P. Benner 4. Stability In contrast to the observations in [12], unstable systems have been encountered when applying BIRKA to industrial problems. Hence, a concept for stability preservation for the reduction with BIRKA is needed. As usual, an (asymptotically) stable linear system will be a system for which it holds Re(λ (A, E)) < 0 for i= 1, … ., n. Therefore, the parametric system in (4) is asymptotically stable if Reðλ ðA þ h N ; EÞÞ < 0 for i= 1, … ., n, and a bilinear system which was i k k k¼1 derived from a parametric system as in (5) is stable if Reðλ ðA þ u N ; EÞÞ < 0 for i k k k¼1 i= 1, … ., n. In general for bilinear systems, the inputs u can be time dependent and the definition of stability for bilinear systems will be stated in terms of the system’s input– output behaviour: Definition 4.1([11]). A bilinear system Ex_ðtÞ¼ AxðtÞþ N u ðtÞxðtÞþ BuðtÞ; k k (11) k¼1 yðtÞ¼ CxðtÞ; is said to be bounded-input-bounded-output (BIBO) stable, if for any bounded input u(t), the output y(t) is bounded on [0, ∞). For the special bilinear systems that result from parametric systems, it is possible to deduce a relation between the eigenvalues of the matrices A and A þ u N . k k k¼1 As N = 0for u resulting from the original linear inputs, only the inputs k k that are time independent will be taken into account and a comparison of the linear and bilinear eigenvalues makes sense. The next Theorem 4.2 and Corollary 4.3 are used to show Proposition 4.4 providing results for the distance between the considered eigenvalues and the stability of the bilinear system in terms of the eigenvalues. nn –1 Theorem 4.2 (Bauer–Fike [20]). If µ is an eigenvalue of M þ S 2 C and X MX = diag(λ , … ., λ ), then 1 n min jλ  μj κ ðXÞjjSjj : (12) i 2 i¼1;...;n –1 nn Corollary 4.3. Let X MX = diag(λ , … ., λ ), and M þ S 2 C . For every eigen- 1 n value λ(M + S) an eigenvalue λ (M) exists such that |λ (M) – λ(M + S)| ≤ κ (X)||S|| . i i 2 2 For a bilinear system with constant inputs u , we derive the following proposition: –1 Proposition 4.4. Let A = X diag(λ , …, λ )X with Re(λ (A)) < – c <0 for all i = 1, … ., 1 n i n. For any j 2 {1, … ., n} there exists an i 2 {1, … ., n} such that jλ ðAÞ λ ðA þ u N Þj<c i j k k k¼1 Mathematical and Computer Modelling of Dynamical Systems 111 if jjujj jjN jj < : (13) 2 2 κ ðXÞ k¼1 In addition Reðλ ðA þ u N ÞÞ <0 for j = 1, …, n. j k k k¼1 Proof. With Corollary 4.3 one concludes: m m X X jλ ðAÞ λ ðA þ u N Þj  κ ðXÞjj u N jj i j k k 2 k k k¼1 k¼1 κ ðXÞjjujj jj N jj 2 k 2 2 k¼1 <c: Assume Reðλ ðA þ u N ÞÞ ≥ 0 for one fixed j 2 {1, … ., n}. As c < |Re(λ (A))| j k k i k¼1 for all i =1, …,n and for j there exists i such that jλ ðAÞ λ ðA þ u N Þj<c one i j k k k¼1 calculates: c< jReðλ ðAÞÞj jReðλ ðAÞÞj þ Reðλ ðA þ u N ÞÞ i j k k k¼1 ¼jReðλ ðAÞÞ  Reðλ ðA þ u N ÞÞj i j k k k¼1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X m 2 m 2 ðReðλ ðAÞ Reðλ ðA þ u N ÞÞ þðImðλ ðAÞÞ  Imðλ ðA þ u N ÞÞ i j k k i j k k k¼1 k¼1 <c; which leads to a contradiction. Therefore, Reðλ ðA þ u N ÞÞ < 0 holds. □ j k k k¼1 For systems with E= I and sufficiently small inputs u and matrices N (cf. (13)), the n k k bilinear system remains stable and every eigenvalue of the bilinear system lies in a neighbour- hood of an eigenvalue of the linear system. For E nonsingular, the Proposition 4.4 remains –1 m valid for E A and u E N . Hence, it will be assumed that the eigenvalues of k k k¼1 m –1 1 1 E A þ u E N and E A are sufficiently close and use stability preserving methods k k k¼1 for linear systems, which works quite well for the considered thermal systems. 4.1. Stability preservation using the systems Gramians For linear systems (i.e. N = 0, k= 1, …, m), stability can be preserved by using the following result due to Yousefi [21]. Basically, Villemagne and Skelton [22] already stated it, whereas Gugercin [23] used it in the context of an interpolatory approach. Yousefi incorporated the fact that the eigenvalues of the reduced model will not exceed a certain value σ. 112 A. Bruns and P. Benner Proposition 4.5 ([21]). Given a linear stable system (A, B, C) with Re(λ (A)) < – σ < 0 for nk T –1 i= 1, …, n. Then for any arbitrary full row rank matrix V 2 R and W = QV(V QV) , T T where Q = Q > 0 satisfies A Q + QA + 2σQ < 0, the reduced model (A,B,C ) is stable r r r and A =W AV satisfies Re(λ (A )) < – σ for i = 1, …, r. r i r For positive semidefinite Q, the Proposition stays valid, if one assumes V QV to be invertible. We generalize it for a system with E ≠ I, E nonsingular, positive semidefinite Q T T and Q ¼ V E QEV nonsingular, which – up to the authors’ knowledge – has not been stated elsewhere. Proposition 4.6. Given a linear stable system (E, A, B, C) with E nonsingular and Re(λ (A, E)) < – σ < 0 for i =1, …, n. Let Q = Q ≥ 0 satisfy T T T A QE þ E QA þ 2σE QE  0: (14) nk T T Then for any arbitrary full row rank matrix V 2 R with Q ¼ V E QEV nonsingular ,A,B,C ) generated with (and therefore Q>0) the reduced model (E r r r r T T 1 W ¼ QEVðV E QEVÞ is stable and satisfies Re(λ (A,E )) ≤ – σ for i =1, … ., r. i r r The proof of the Proposition follows exactly the proof of Yousefi (cf. Proposition 4.5). However, as we have introduced two generalizations – the presence of the E matrix and the nonstrict Lyapunov inequality (cf. Equation (14)) – we state it here for completeness. Proof. Multiplying Equation (14) with V and V and making use of T T 1 T T T T T T T T I ¼ðV E QEVÞ ðV E QEVÞ ¼ðV E QEVÞ ðV E QEVÞ; leads to: T T T T T T V A QEV þ V E QAV þ 2σV E QEV  0 T T 1 T T T T T T T T T T ) V A QEVðV E QEVÞ ðV E QEVÞ ðV E QEVÞ ðV E QEVÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ^ W T 1 T T T T T T T T T T þ V E QEVðV E QEVÞ V E QEV ðV E QEVÞ V E Q AV |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} W ^ W T T 1 T T T T T T T T T T þ 2σV E QEVðV E QEVÞ ðV E QEVÞ ðV E QEVÞ V E Q EV  0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} W ^ W T T T T T T T T T ^ ^ ^ ) V A WQW EV þ V E WQW AV þ 2σV E WQW EV  0 T T T ^ ^ ^ ) A QE þ E QA þ 2σE QE  0 r r r r r r ^ ^ )ðA þ σE Þ QE þ E QðA þ σE Þ 0: r r r r r r Mathematical and Computer Modelling of Dynamical Systems 113 T T T T T r Using the identity E ¼ W EV ¼ðV E QEVÞ V E QEV ¼ I , let λ and v be any r r i eigenvalue and eigenvector of [A + σI ], then: r r T T ^ ^ ^ ^ ðA þ σI Þ Q þ QðA þ σI Þ 0 ) v ðA þ σI Þ Qv þ v QðA þ σI Þv  0 r r r r r r i r r i i i r r ^ ^ ) λ v Qv þ λ v Qv  0 i i i i i i r r )ðλ þ λ Þv Qv  0 i i i ) 2Reðλ Þv Qv  0 i i ðv Qv >0Þ) Reðλ Þ 0: i i Consider the eigenvalues of the reduced system (which are the eigenvalues of A as r r r E = I ) With λ v ¼ðA þ σI Þv this leads to A v ¼ λ v  σv ¼ðλ  σÞv.As r r i r r i r i i i i i i i r r Reðλ Þ 0 and – σ < 0, one can conclude that Reðλ  σÞ < 0 and therefore that the i i reduced system is stable. □ The dual result is also true: Proposition 4.7. Given a linear stable system (E, A, B, C) with E nonsingular and Re(λ (A, E)) < – σ < 0 for i = 1, …, n. Then for any arbitrary full row rank matrix nr T W 2 R and P = P ≥ 0 which satisfies T T T APE þ EPA þ 2σEPE  0 (15) T T and nonsingular P ¼ WEPE W , the reduced model (E,A,B,C ) generated with r r r r T T T T V ¼ PE W ðWEPE W Þ is stable and satisfies Re(λ (A,E )) < – σ for i = 1, …, r. i r r Proof. The proof is analogous to the one of Proposition 4.6. □ For the calculation of the projection matrix W the following Lyapunov equation is solved: T T ðA þ σEÞ QE þ E QðA þ σEÞ¼C C  0; (16) for a σ < |Re(λ (A, E))|. This leads to the calculation of the projection matrix max T T W ¼ QEVðV E QEVÞ as in Proposition 4.6. The solution of the Lyapunov Equation (16) is positive semidefinite, as the shifted system (A + σE, E) remains asymptotically stable. Equation (16) can be solved using the low rank alternating direction implicit (ADI) iteration (cf. for example [24, 25]) which generates a low rank factor Z, such that Q ≈ Z Z. The calculated low rank 114 A. Bruns and P. Benner T T T matrix Q  V E Z ZEV can be singular. This always occurs in the case where rk(Z) <rk T T 2 (V) =r. Even if rk(V) ≥ rk(Z) one cannot conclude that V E QEV is nonsingular , but for rk(V) relatively small compared to rk(Z) it is often the case. The proposed stability preservation method will be used for one of our models in Section 7.1.1. Solving large Lyapunov equations is numerically demanding. For large systems (n > 500, 000) it might be impossible – even with highly developed methods such as the ADI algorithm. Hence, this stability preserving method will reach its limitations when the systems dimensions are getting too large. 4.2. Stability preservation via one-sided projections In the special case of symmetric matrices E, A and N and positive definite E, another possibility for preserving stability is using only a single projection matrix (one-sided method). The matrices of the thermal system provided in Equation (24) in Section 6 have exactly these properties, and therefore this stability preservation is of interest. Proposition 4.8 ([26]). Given a linear system (i.e. N =0) with A, E symmetric, the system T T is asymptotically stable if and only if E = E >0 and A = A <0. nn nr Corollary 4.9 (4.4 [18],). Let A 2 R be a symmetric matrix. Let V 2 R have orthonormal columns. Let A = V AV. Then λ ðAÞ λ ðA Þ λ ðAÞ; i ¼ i; ... ; r (17) i i r iþnr T T Corollary 4.10 ([26]). Given a linear system with E = E >0 and A = A with λ (A) < 0, nr T T for i = 1, …, n. Let V 2 R have orthonormal columns. Let A =V AV and E =V EV. r r Then the reduced system is asymptotically stable. Hence for linear systems with A, E symmetric and E positive definite stability can be preserved via one-sided projections. Changing to bilinear systems, one needs the follow- ing results. nn Proposition 4.11 (Weyl [18], ). Let A; B 2 R be two symmetric matrices. Let λ (A) and λ (B) for i = 1, …, n be the eigenvalues of A and B with λ ðAÞ    λ ðAÞ and i 1 n λ ðBÞ    λ ðBÞ. Then it holds: 1 n λ ðA þ BÞ2½λ ðAÞþ λ ðBÞ; λ ðAÞþ λ ðBÞ for i ¼ 1; ... ; n (18) i i n i 1 Corollary 4.12 ([18]). Under the assumptions of Proposition 4.11 it holds jλ ðA þ BÞ λ ðAÞj  jjBjj for i ¼ 1; ... ; n (19) i i As in the context of the Bauer–Fike Theorem 4.2, the eigenvalues of a bilinear system, derived from a linear, parametric system can be related to the eigenvalues of the linear system. Using the previous results, we obtain the following Corollaries: Mathematical and Computer Modelling of Dynamical Systems 115 nn nn Corollary 4.13. Let u 2 R for k = 1, …, m, A 2 R and N 2 R symmetric with k k nr eigenvalues 0> λ ðAÞ    λ ðAÞ and λ ðN Þ λ ðN Þ. Given that V 2 R 1 n 1 k n k T T has orthonormal columns and A =V AV and (N ) =V N V, then it holds r r k k m m X X jλ ðA þ u ðN Þ Þ λ ðA Þj  jjujj jjN jj : (20) i r k r i r k k 2 2 k¼1 k¼1 Proof. Corollary 4.12 leads to m m X X jλ ðA þ u ðN Þ Þ λ ðA Þj  jj u ðN Þ jj i r k r i r k r k k 2 k¼1 k¼1 jjujj jjðN Þ jj : 2 k 2 k¼1 As N and (N ) are symmetric, they are normal and therefore fulfil ||(N ) || = max k r k r k 2 i =l, |λ ((N ) ) | = max{|λ ((N ) )|,|λ ((N ) )|} and ||N || = max |λ (N )|= max …, r i r k 1 r k r r k k 2 i =1, …,r i k {|λ (N )|, |λ (N )|}. With Corollary 4.9 one concludes λ ðN Þ λ ððN Þ Þ λ ðN Þ. This 1 k i r n k 1 k r k leads to jjðN Þ jj jjN jj and therefore Equation (20) holds. □ r k k 2 2 In addition one obtains: Corollary 4.14. Under the assumptions of Corollary 4.13 let c 2 R with c < P P m m jλ ðA Þj ¼ jλ ðA Þj. If jjujj jjN jj <c then λ ðA þ u ðN Þ Þ < 0. max r 1 r k i r k r 2 k¼1 2 k¼1 k Proof. Assume λ ðA þ u ðN Þ Þ 0 and calculate using Equation (20): i r k r k¼1 k c< jλ ðA Þj<jλ ðA Þj 1 r i r jλ ðA Þj þ λ ðA þ u ðN Þ Þ i r i r k r k¼1 ¼jλ ðA Þ λ ðA þ u ðN Þ Þj <c: i r i r k r k¼1 This leads to a contradiction and therefore λ ðA þ u ðN Þ Þ < 0. □ i r k r k¼1 k If jjujj jjN jj is sufficiently small, one can assume that the eigenvalues 2 2 k¼1 P P m m λ ðA þ u ðN Þ Þ and λ (A ) are similar. In addition, λ ðA þ u ðN Þ Þ < 0. i r k r i r i r k r k¼1 k k¼1 k Hence, starting from a model with λ (A) < 0 for all i =1, …, n and jjujj jjN jj i k 2 k¼1 2 sufficiently small leads to stable reduced bilinear systems. 4.3. Stability preservation – the workflow As the reduced models, that have been calculated with the stabilization process using the Gramians, are in most cases better than those generated by a one-sided approach, the workflow in Figure 1 applies. The reader should note that l r indicates the fact, that the matrix Q (cf. Section 4.1) can always be singular, but for the case of l r, it is more likely that Q is invertible. 116 A. Bruns and P. Benner Fix a reduced order r. Is solving equation (16) by an ADI iteration possible? Yes No Solve the equation (16) and deter- Use a one-sided approach. mine the rank rk(Z)= l. l  r l < r Stability preservation via Proposition 4.6. Figure 1. Proposed workflow for stabilization. 5. Heat transfer modelling 5.1. Mode of operation of an electrical motor An electrical motor converts electrical energy into mechanical work, which is produced by the interaction of an electrical current and a magnetic field. One part of the motor – the so- called stator – consists of several coils wound around an iron core. When a voltage is applied, a current is induced in the coil. Inside the counterpart – the so-called rotor – a magnetic field is generated either by a permanent magnet or by an electromagnet. The interaction of this magnetic field with the current in the stator results in a rotation of the rotor. See Figures 2 and 3. Actuating the motor with electrical currents leads to an increase in temperature in its different components due to thermal losses. It is important to analyse the influence of this temperature change on the materials of the motor, as it affects its lifespan. This is achieved by carrying out a thermal analysis. 5.2. Thermal modelling of an electrical motor For a thermal analysis, several physical effects have to be considered and can be modelled based on the three main types of heat transfer: heat conductance, convection and radiation. magnets coil stator rotor Figure 2. Drawing of a slice through an electrical motor. Mathematical and Computer Modelling of Dynamical Systems 117 (a) Drive unit and generator in one: the (b) Generator for commercial vehicles. The Bosch integrated motor generator. configuration of the coils is the same as for an electrical motor Figure 3. Two Bosch components illustrating the structure of an electrical motor. Pictures: Bosch. © [Bosch]. Reproduced by permission of Robert Bosch GmbH. For a broad overview of heat and mass transfer, see for example the book of Baehr and Stephan [27]. The main heat sources in the electrical motor are thermal losses, resulting from the current in the coil of the stator and/or rotor. The motor has to fulfil various operational requirements and therefore different current profiles have to be considered. The tempera- ture on certain parts of the motor (for example, the flange) should not exceed a specified upper limit, because these parts are in contact with other temperature sensitive compo- nents. This upper limit is built into the model as a fixed temperature. The motor is surrounded by air, therefore convection has to be considered. In general, convection describes the heat transfer between a fluid (here: surrounding air) and a solid. This effect can be modelled via a so-called heat transfer coefficient h. Different environ- ments have to be analysed, as the motor needs to work in a large temperature range (arctic winter, tropical summer); therefore, different ambient temperatures are examined in the model. Varying the heat transfer coefficients represent different cooling strategies or different interaction scenarios of the motor with its environment. Various parts of the motor are not directly attached and the resulting thermal resistance has to be modelled by a contact heat transfer coefficient. Varying this parameter, the small gap between the two materials can be considered as filled with air or an insulation material. The motor is built from various materials such as steel, copper and plastics. These materials have different properties, among others the density ρ, the specific heat C and the thermal conductivity k. Here, these material parameters will not be varied. As the motor temperature remains relatively small, the effect of radiation is not of great impor- tance, and therefore it will not be analysed. 6. Model parametrization To determine the thermal behaviour of the motor while in operation, the temperature field T(x, t) (dependent on location x and time t) has to be examined. The temperature field T(x, t) within a domain Ω 2 R for times t 2½ 0; t can be calculated by using the heat end equation @TxðÞ ; t ρC ¼ kTxðÞ ; t þSxðÞ ; t : (21) @t 118 A. Bruns and P. Benner using constant material properties ρ,C,k and a heat source S, which in case of an electrical motor corresponds to the thermal losses of the coil of stator or rotor. The derivation of the heat equation can be found, e.g. in [27]. On interfaces and outer surfaces, now called boundaries and denoted as Γ 2 R , different conditions have to be given, depending on the situation of interest. They are mathematically formulated as follows: ● Dirichlet boundary conditions: Tðx; tÞ¼ T ðtÞ on the boundary Γ : D D These conditions correspond to fixed temperatures on surfaces. Neumann boundary conditions: @TxðÞ ; t k ¼ q_ in Γ ; N N @n where q_ is a given heat flux on the boundary. In the considered model of the electrical motor, this kind of boundary condition does not appear. Robin boundary conditions: @TxðÞ ; t k hTðÞ  T in Γ ; 1 R @n where h denotes the heat transfer coefficient. The temperature T is the temperature of the surrounding domain. Interface conditions: A thermal resistance between two surfaces can be modelled on the interface by a thermal contact conductance coefficient h . The interface I will be considered as two surfaces: I with temperature T and I with temperature T . The following equation 1 1 2 2 applies: @TjIðÞ x; t @TjIðÞ x; t 1 2 k ¼k ¼ h TxðÞ ; tj TxðÞ ; tj : 2 1 c I I 1 2 @n @n A spatial discretization of the heat Equation (21) is then performed by using the finite element method. Incorporating the boundary and interface conditions and using a finite element basis function ψ; ψj ¼ 0, the weak formulation reads as follows: Z Z Z @TxðÞ ; t ψðÞ x ρC dx þ ψðÞ x  kTxðÞ ; t dx þ ψðÞ x hTðÞ x; t ds @t Ω Ω Γ Z Z þ ψðÞ x h TxðÞ ; t ds  ψðÞ x h TxðÞ ; t ds (22) c c I I 1 2 Z Z Z ¼ ψðÞ x SxðÞ ; t dx þ ψðÞ xðÞ qN _ ds þ ψðÞ x hT ds: Ω Γ Γ N R The material parameters ρ,C and k are considered as constant. With finite element basis functions ψ ðxÞ the temperature is approximated as follows, k Mathematical and Computer Modelling of Dynamical Systems 119 Tðx; tÞ T ðtÞψ ðxÞ: k¼1 This leads to the following discretized equation: 2 3 SðtÞ hT 6 7 ETðtÞ¼ ðA þ hN þ h N ÞTðtÞþ B: : (23) 1 c 2 4 5 q_ R R R R As ψ ðxÞψ ðxÞdx ¼ ψ ðxÞψ ðxÞdx and ψ ðxÞ ψ ðxÞdx ¼ ψ ðxÞ ψ ðxÞdx; k j j k k j j k the matrices E, A and N for the considered class of systems are symmetric and E is in addition positive definite. One should note that a thermal conductance can lead to a singular matrix A, whereas A + h N is nonsingular. Hence a new parameter h ¼ h  s is introduced, and c 2 c c ~ ~ A þ h N with A ¼ A þ sN andafixedshift s will be considered. c 2 2 7. Results The thermal analysis is carried out using COMSOL Multiphysics, version 3.5a [28]. By exporting several matrices from COMSOL and a thorough analysis of the underlying equations, it is possible to reconstruct a parametric model with variable parameters and loads. 2 3 h T 1 1 6 7 ! 6 7 q . v 6 7 X X 6 7 _ ~ ~ ETðtÞ¼ A þ h N þ ðh Þ N TðtÞþ B  ; 6 7 i i c k k h T q 1 6 7 (24) i¼1 k¼qþ1 6 7 4 5 LðtÞ yðtÞ¼ CTðtÞ; where q is the number of heat transfer coefficients h,and v – q is the number of contact heat transfer coefficients h .In Figure 4, the modelled motor part is shown. One can see parts of stator, coil, housing and some insulation parts. The following loads and parameters need to be considered. On top of the housing a temperature T is fixed to include a specified maximum temperature. The coil losses L(t) are incorporated into the coil. Heat transfer by convection is modelled at seven different locations, for example on coil and housing, resulting in 7 heat transfer coefficients (i.e. q = 7). Thermal resistance is incorporated at six different locations, for example between the insulation and the stator or the insulation and the coil (i.e. v= 13). The size of the model is n = 41, 199 and the original transient COMSOL simulation for one parameter setting takes about 90 minutes. There are two different models of the electrical motor that have been examined. The first one is simulated using only heat transfer coefficients as parameters and ignoring the effects of thermal resistance between some motor parts. This leads to a model with 7 120 A. Bruns and P. Benner Figure 4. The Comsol model simulates the heat transfer in a stator slice, without the rotor. parameters and 4 loads. The second model simulates the thermal resistances and therefore contact heat transfer coefficients have to be considered. Hence, the resulting model has 13 parameters and 4 loads. Each of the resulting parametric systems (24) is reformulated as a bilinear system by using the procedure explained in [7] and afterwards reduced using BIRKA (Algorithm 1) with the calculation of the projection matrices V and W as explained in (9). The infinite sum is truncated after 10 summands. The calculations were performed using MATLAB [19] on 12CPUs with 3GB RAM each. 7.1. Model 1 – no contact heat transfer coefficients 7.1.1. General results The first reduction has been executed using the stability preservation methods in Section 4.1. Hence, the stability of the original model is preserved by calculating the projection matrix W as proposed in Proposition 4.6. The reduction took about 11 hours. This –7 corresponds to 16 iterations and a change in the eigenvalues of less than 10 . The reduced model has order r = 50 and can be simulated within 10 seconds, which corresponds to a speed-up of over 500 compared to the original simulation time of 90 minutes. Compared to the reduction time of 11 hours, the original model could have been simulated about 8 times. As an optimal parameter set for predefined operating conditions is determined within an optimization process, it requires up to 10, 000 simulations of the model – the reduced model pays off. Now the results at four different locations will be examined. The temperatures at the front of the stator, at the coil and at two different points on the insulation have to be monitored. The original solution is compared to the solution which was computed using the reduced model. Figure 5 shows that there is only a small deviation which can be seen in the error plots. The relative error in temperature is smaller than 0.07 K, corresponding to a relative maximal error of less than 2  10 . It is important to make sure that the reduced model gives reliable results over a wide range of parameter values and inputs. The heat transfer coefficients in our model range from 5 to 100, and the coil losses L(t) and the ambient temperature T also need to be varied. For all these variations the reduced ∞ Mathematical and Computer Modelling of Dynamical Systems 121 Figure 5. Different outputs – original model compared to reduced model and relative and absolute errors. model gives an excellent approximation of the full model. Take for example Figure 6, where the behaviour of the temperature for six different heat transfer coefficients on the coil is shown. The error plots show that again the relative and absolute errors are sufficiently small. 122 A. Bruns and P. Benner Figure 6. Temperature curve for six different values (5, 25, 45, 65, 85, 100[W/m K]) of the heat transfer coefficient on the coil, and the relative and absolute errors. 7.1.2. Stability preserving – comparing the different approaches As explained in Sections 4.1 and 4.2, stability can be preserved by different procedures. Here, the following approaches will be examined: gramian BIRKA: The reduced model is calculated by using V as in Algorithm 1 with Equation (9), and the matrix W is calculated using Theorem 4.6. Mathematical and Computer Modelling of Dynamical Systems 123 BIRKA-tS: The projection matrices V and W are calculated with Algorithm 1. Stability is not preserved. Hence, in every step of the iterative process the generated reduced system is saved, and a stable system is chosen from within these systems. This stable system does not necessarily exist, and even if, it is possibly not an optimal reduced system, as it is not necessarily the final reduced system. BIRKA-oS: The reduced model is calculated with Algorithm 1. Only in the last step a one- sided projection with V is used. only V: Calculate V as in Algorithm 1. In every step of the algorithm a one-sided projection is used to calculate the reduced model. For the outputs on top of the coil and on the insulation between coil and stator, all reduced models of order 40–60 have a sufficient accuracy. Figure 7 shows the relative and absolute errors of the outputs 1 and 4 on the stator front and the insulation on the top of the stator, respectively. For the reduced order 40, the original Algorithm 1 (cf. dash-dotted line in Figure 7) performs well, whereas the other approaches are still not satisfying. However, if the order is increased, the approximation improves. With the exception of the approach of only using V within the algorithm. In Figure 8, the approaches BIRKA-oS and only V are compared for different orders. For the latter one finds, that the order has to be increased up to r = 100 to get results, that are as good as for the other methods. Figure 7. Output 1 and 4 for different reduced orders error plots. 124 A. Bruns and P. Benner Figure 8. One-sided methods. 7.2. Model 2 – contact heat transfer coefficients In the second model, thermal resistance has been taken into account. Six additional parameters (contact heat transfer coefficients h ) are incorporated into the model. They W W range from 200 up to 3600 . As mentioned in Section 6, these parameters can lead 2 2 m K m K to a singular matrix A, and a shift s needs to be introduced leading to a nonsingular matrix min max A ¼ A þ sN. For every given h 2½h ; h , one shift (the centre of the interval) is c c needed for achieving a nonsingular A. Making use of the shifts might affect the stability of Mathematical and Computer Modelling of Dynamical Systems 125 ~ ~ ~ the reduced model. For A stable, one cannot guarantee that A  sN ¼ A remains stable. r r r r The stability of A and A can be connected using the statements derived in Section 4.2 and the Bauer–Fike Theorem 4.2. For this model, the stability preservation using Theorem 4.6 will fail, because the size of a reduced model will be larger than the rank of the low rank factor in the ADI iteration. Hence, stability can only be preserved using a one-sided projection. This leads to larger reduced orders than an original BIRKA would lead to. For the reduction, the following approaches are used: ● only V: The projection matrix V is calculated as in Algorithm 1. In every step of the algorithm a one-sided projection is used to calculate the reduced model. ● BIRKA-tS: The projection matrices V and W are calculated with Algorithm 1 out of a reduced model generated by only V. Stability is not preserved. Hence, in every step of the iterative process the generated reduced system is saved, and a stable system is chosen from within these systems. This stable system does not necessarily exist and, even if, it is possibly not an optimal reduced system, as it is not necessarily the final reduced system. BIRKA-oS: The reduced model is calculated with Algorithm 1 out of a reduced model generated by only V. Only in the last step a one-sided projection with V is used. The reduction was performed using the one-sided approach only V and took about 3 days and 3 hours. The reduced model has order r = 600 and can be simulated within 60 seconds, which corresponds to a speed-up of about 90 compared to the original simulation time of 90 minutes. This reduced model leads to a good approximation of the original model within the whole parameter range. This can be seen for instance in Figure 9, where the variation of the heat transfer coefficient on the coil is shown. The two approaches BIRKA-tS and BIRKA-oS use the reduced model calculated with only V and reduce it again. This leads to models of order r = 300, and additional 12 hours for the reduction process. These models can be simulated in 15 seconds, which corresponds to a speed-up of over 300 compared to the original simulation time. Figure 10 shows the original and the reduced models generated by using the above approaches, and the largest error which occurs in output 3. The reduced models generated with only V and BIRKA- tS show sufficient accuracy, whereas the approximation with BIRKA-oS is still not accurate enough. 8. Discussion and conclusion While applying BIRKA to thermal industrial models, several problems were encoun- tered. First, approximating the original calculation of the projection matrices as shown in (9) could lead to a divergence of the algorithm. This divergence can be detected by the procedure shown in Section 3; however, there still might be other causes for divergence. Second, one class of physical parameters (contact heat transfer coefficients) leads to a singular stiffness matrix A. Hence, the reduction is performed using a shifted matrix A, as explained in Section 6. Third, stability is not preserved by the algorithm and two main strategies to stabilize the reduced models are presented in Sections 4.1 and 4.2. Considering the speed-up and the accuracy of the reduced model, BIRKA seems to be an option for the reduction of industrial thermal models. As the reduction times are relatively long (especially for the second model, cf. Section 7.2), the reduction is only reasonable if a large amount of test cases with different parameter samples need to be simulated. Other methods for the reduction of parametric models use a different approach. In general, models are reduced in several parameter points, and an 126 A. Bruns and P. Benner Figure 9. Temperature curve for six different values (5, 25, 45, 65, 85, 100[W/m2K]) of the heat transfer coefficient on the coil, and the relative and absolute errors. interpolation procedure is used to obtain reduced models in new parameter points (as for example in [2–6, 29]). For models with a small parameter range, these methods might lead to smaller reduced orders, as BIRKA offers a good accuracy in a large parameter space, which might not be needed for the considered model. In addition, BIRKA is restricted to the special parameter dependence given in Equation (4). Linear models with a parameter dependence in the matrices E, B or C or with a nonaffine parameter dependence cannot be reduced directly with this algorithm, and hence other methods need to be applied. The methods given in [2–6, 29]are basedon linear Mathematical and Computer Modelling of Dynamical Systems 127 Figure 10. Reduction of model 2 with different approaches, largest error – output 3. reduction in certain parameter points. If the reduction in one sample point is not time- consuming and a small number of sampling points is needed, it is likely that the method outperforms BIRKA in terms of reduction time. However, these methods can only guarantee a small reduction error in the neighbourhood of the sampling points. BIRKA, in contrast, will yield a model with a globally small reduction error. 128 A. Bruns and P. Benner Notes 1.kk M # M ¼kk M kk M ½ 30 , Corollary 13.11. 1 2 1 2 2 2 2 0 1 0 1 0 1 1 100 00 B C T T T @ AB C @ A 2. Take n ¼ 4; E ¼ I ; V E Z ¼ 0 011 ¼ 00 @ A 0 1 1 000 10 References [1] A.C. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, Philadelphia, 2005. [2] D.S. Weile, E. Michielssen, E. Grimme, and K. Gallivan, A method for generating rational interpolant reduced order models of two-parameter linear systems, Appl. Math. Letters. 12 (1999), pp. 93–102. doi:10.1016/S0893-9659(99)00063-4 [3] L. Daniel, O. Siong, L.S. Chay, K.H. Lee, J. White, and J. White, A multiparameter moment matching model reduction approach for generating geometrically parameterized interconnect performance models, IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 23 (2004), pp. 678–693. doi:10.1109/TCAD.2004.826583 [4] U. Baur, C.A. Beattie, P. Benner, and S. Gugercin, Interpolatory projection methods for parameterized model reduction, SIAM J. Sci. Comp. 33 (2011), pp. 2489–2518. doi:10.1137/090776925 [5] H. Panzer, J. Mohring, R. Eid, and B. Lohmann, Parametric model order reduction by matrix interpolation,at – Automatisierungstechnik. 58 (2010), pp. 475–484. doi:10.1524/ auto.2010.0863 [6] D. Amsallem and C. Farhat, An online method for interpolating linear parametric reduced- order models, SIAM J. Sci. Comp. 33 (2011), pp. 2169–2198. doi:10.1137/100813051 [7] P. Benner and T. Breiten, On H -model reduction of linear parameter-varying systems,in Proceedings of Applied Mathematics and Mechanics, WILEY-VCH Verlag, Weinheim, 2011, pp. 805–806. doi:http://dx.doi.org/10.1002/pamm.201110391 [8] S. Al-Baiyat and M. Bettayeb, A New Model Reduction Scheme for k-Power Bilinear Systems, Proceedings of the 32nd IEEE Conference on Decision and Control Vol. 1, 1993, pp. 22–27. [9] Z. Bai and D. Skoogh, A projection method for model reduction of bilinear dynamical systems, Linear Algebra Appl. 415 (2006), pp. 406–425. doi:10.1016/j.laa.2005.04.032 [10] T. Breiten and T. Damm, Krylov subspace methods for model order reduction of bilinear control systems, Syst. Control Lett. 59 (2010), pp. 443–450. doi:10.1016/j. sysconle.2010.06.003 [11] L. Zhang and J. Lam, On H model reduction of bilinear systems, Automatica. 38 (2002), pp. 205–216. doi:10.1016/S0005-1098(01)00204-7 [12] P. Benner and T. Breiten, Interpolation-based H -model reduction of bilinear control systems, SIAM J. Matrix Anal. Appl. 33 (2012), pp. 859–885. doi:10.1137/110836742 [13] Y. Bertin, E. Videcoq, S. Thieblin, and D. Petit, Thermal behavior of an electrical motor through a reduced model, IEEE T. Energy Conver. 15 (2000), pp. 129–134. doi:10.1109/ 60.866989 [14] Z. Gao, R. Colby, T. Habetler, and R. Harley, A model reduction perspective on thermal models for induction machine overload relays, IEEE Trans. Ind. Electronics. 55 (2008), pp. 3525–3534. doi:10.1109/TIE.2008.926772 [15] K. Zhou, J. Pries, H. Hofmann, Y. Kim, T.-K. Lee, and Z. Filipi, Computationally-efficient finite-element-based thermal models of electric machines, Vehicle Power and Propulsion Conference (VPPC), 2011 IEEE, Chicago, IL, IEEE, 6–9 September 2011, 1–6. Available at http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6043205&isnumber=6042961 (Accessed 2 June 2014). [16] G.M. Flagg, Interpolation Methods for the Model Reduction of Bilinear Systems, Ph.D. thesis, Virginia Polytechnic Institute and State University, Blacksburg, VI, 2012. [17] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2003. [18] G.W. Stewart and J.G. Sun, Matrix Perturbation Theory (Computer Science and Scientific Computing), Academic Press, San Diego, CA, 1990. [19] MATLAB, version 7.10.0 (R2010a), The MathWorks Inc., Natick, Massachusetts, 2010, ©The MathWorks, Inc. MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See www.mathworks.com/trademarks for a list of additional trademarks. Mathematical and Computer Modelling of Dynamical Systems 129 Other product or brand names may be trademarks or registered trademarks of their respective holders. [20] G. Golub and C.V. Loan, Matrix Computations, 3rd ed. John Hopkins University Press, Baltimore, MD, 1996. [21] A. Yousefi, Preserving stability in model and controller reduction with application to embedded systems, Ph.D. thesis, Technische Universität München, München, 2006. [22] C.D. Villemagne and R.E. Skelton, Model reductions using a projection formulation,J. Control. 46 (1987), pp. 2141–2169. doi:10.1080/00207178708934040. [23] S. Gugercin, An iterative SVD-Krylov based method for model reduction of large-scale dynamical systems, Linear Algebra Appl. 428 (2008), pp. 1964–1986. doi:10.1016/j.laa.2007.10.041. [24] P. Benner, J.R. Li, and T. Penzl, Numerical solution of large-scale lyapunov equations, riccati equations, and linear-quadratic optimal control problems, Numer. Linear Algebra Appl. 15 (2008), pp. 755–777. doi:10.1002/nla.622. [25] J. Saak, Efficient numerical solution of large scale algebraic matrix equations in PDE control and model order reduction, Ph.D. thesis, Technische Universität, Chemnitz, 2009. [26] R. Castan Selga, The matrix measure framework for projection-based model order reduction, Ph.D. thesis, Technische Universität, München, 2011. [27] H.D. Baehr and K. Stephan, Heat and Mass Transfer, Springer, Berlin, 2006. [28] COMSOL-Multiphysics, version 3.5a, COMSOL is a registered trademark of COMSOL AB., [29] P. Benner and L. Feng, A robust algorithm for parametric model order reduction based on implicit moment matching, Reduced Order Methods for Modeling and Computational Reduction, in MS&A, A. Quarteroni and G. Rozza, eds., Springer, 2014, vol. 9, 159–186. [30] A.J. Laub, Matrix Analysis for Scientists and Engineers, SIAM, Philadelphia, 2005. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

Parametric model order reduction of thermal models using the bilinear interpolatory rational Krylov algorithm

Download PDF
27 pages

Loading next page...
 
/lp/taylor-francis/parametric-model-order-reduction-of-thermal-models-using-the-bilinear-iINJQqeJBP

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Taylor & Francis
Copyright
© 2014 Taylor & Francis
ISSN
1744-5051
eISSN
1387-3954
DOI
10.1080/13873954.2014.924534
Publisher site
See Article on Publisher Site

Abstract

Mathematical and Computer Modelling of Dynamical Systems, 2015 Vol. 21, No. 2, 103–129, http://dx.doi.org/10.1080/13873954.2014.924534 Parametric model order reduction of thermal models using the bilinear interpolatory rational Krylov algorithm a b Angelika Bruns * and Peter Benner a b Robert Bosch GmbH, CR/ARH2, Robert-Bosch-Platz 1, 70839 Gerlingen, Germany; Max Planck Institute for Dynamics of Complex Technical Systems, Computational Methods in Systems and Control Theory, Sandtorstr. 1, 39106 Magdeburg, Germany (Received 25 September 2013; accepted 12 May 2014) The Bilinear Interpolatory Rational Krylov Algorithm (BIRKA; P. Benner and T. Breiten, Interpolation-based H -model reduction of bilinear control systems, SIAM J. Matrix Anal. Appl. 33 (2012), pp. 859–885. doi:10.1137/110836742) is a recently developed method for Model Order Reduction (MOR) of bilinear systems. Here, it is used and further developed for a certain class of parametric systems. As BIRKA does not preserve stability, two different approaches generating stable reduced models are pre- sented. In addition, the convergence for a modified version of BIRKA for large systems is analysed and a method for detecting divergence possibly resulting from this modifica- tion is proposed. The behaviour of the algorithm is analysed using a finite element model for the thermal analysis of an electrical motor. The reduction of two different motor models, incorporating seven and thirteen different physical parameters, is performed. Keywords: parametric model order reduction; stability preservation; thermal heat transfer; finite element modelling Subject Classification: 34K17; 93A15 1. Introduction The design of a new product in industry is often accompanied by computer simulations, which are an excellent tool for analysing the behaviour of the involved components. For the estimation of the product’s lifespan, thermal simulations are employed to detect whether an excessive increase in temperature will lead to a fatigue or even deterioration of the materials. Thermal simulations are often carried out using the finite element method, for which excellent tools exist. After a finite element discretization, the tool solves a system of differential equations of the following form: Ext _ðÞ ¼ AxðÞ t þ BuðÞ t ; (1) ytðÞ ¼ CxðÞ t ; n nn nm pn with the state vector x(t) 2 R , the constant matrices E; A 2 R , B 2 R , C 2 R , m p the inputs uðtÞ2 R , and outputs yðtÞ2 R . Usually, the thermal finite element models are complex (often more than 500,000 degrees of freedom), and require therefore a large numerical effort to converge to a solution. *Corresponding author. Email: angelika.bruns@de.bosch.com © 2014 Taylor & Francis 104 A. Bruns and P. Benner Model Order Reduction (MOR) is a powerful method to reduce the dimension of the underlying system and therefore the simulation time significantly. First, an approximation nr of the original state vector xtðÞ  Vx ðÞ t with V 2 R is required while the dimension r of x (t) should be much smaller than the original dimension n. Equation (1) now reduces to EVx_ ðÞ t ¼ AVx ðÞ t þ BuðÞ t þ "ðÞ t ; r r (2) y ðÞ t ¼ CVx ðÞ t ; r r with a residual ε(t). By multiplication from the left with the transpose of a matrix nr W 2 R which is chosen such that the residual ε(t) is orthogonal to the subspace spanned by its columns (i.e. W ε(t) = 0), the following reduced system is obtained: E A B r r r zfflfflffl}|fflfflffl{ zfflfflffl}|fflfflffl{ zffl}|ffl{ T T T W EV x_ ðÞ t ¼ W AV x ðÞ t þ W B utðÞ; r r (3) y ðÞ t ¼ CV x ðÞ t ; r r |{z} rr rm pr p1 where E , A 2 R , B 2 R , C 2 R and y ðtÞ2 R . The goal of projection- r r r r r based MOR is to determine the matrices V and W. For this class of systems, there exist many well-known methods such as Balanced Truncation, Proper Orthogonal Decomposition (POD) or Krylov subspace methods. An overview of these methods can be found in [1]. In the context of thermal simulations, several physical effects have to be considered. Heat sources, convection between the component and the surrounding environment and predefined temperatures on certain parts have to be taken into account. For this reason, the system (1) has to be reformulated to include the parameters h and the loads u(t) corresponding to the considered physical effects above: Ext _ðÞ ¼ A þ h N xtðÞþ BuðÞ t ; i i (4) i¼1 ytðÞ ¼ CxðÞ t ; nn nl kn n l k with E, A, N 2 R , B 2 R , C 2 R , xðtÞ2 R , uðtÞ2 R and yðtÞ2 R . Different parameters and loads denote different situations and designs. It is therefore desirable to reduce the system while keeping the dependence on the parameters h . This problem belongs to the class of parametric systems, for which reduction – the so-called parametric Model Order Reduction (pMOR) – was examined by a number of authors. Using a Moment-Matching approach, Weile et al. [2] reduced systems with two parameters, while Daniel et al. [3] extended this method to the multiple parameter case. Another well-known idea is the calculation of different reduced models for different parameter values followed by an interpolation between their underlying matrices. For example, Baur et al. [4] combine this idea with interpolatory MOR, Panzer et al. [5] use Krylov subspace methods and Amsallem et al. develop an interpolation procedure on a certain mani- fold [6]. Recently, Breiten and the second author proposed to rewrite a parametric system (4) as a bilinear system [7] of the form: Mathematical and Computer Modelling of Dynamical Systems 105 Ext _ðÞ ¼ AxðÞ t þ N u xtðÞþ BuðÞ t ; i i i¼1 (5) bil ytðÞ ¼ CxðÞ t ; Instead of reducing a linear, parametric system, the reduction of a bilinear system with parameters as inputs can now be realized. Al-Baiyat and Bettayeb studied the reduction of bilinear systems based on balancing, see for example [8], whereas Bai and Skoogh developed methods for reducing them with Krylov subspace methods [9]. This approach was further examined by Breiten and Damm [10]. An H optimal method was developed by Zhang and Lam [11] and transferred to interpolatory MOR by Breiten and the second author [12]. This new algorithm, the so-called Bilinear Interpolatory Rational Krylov Algorithm (BIRKA), will be examined within an industrial context in this work and will be applied to the reduction of an electrical motor. Other methods for the efficient simulation of electrical motor models can be found in the following works: Bertin et al. [13] and Gao et al. [14] apply their reduction methods to thermal networks. As it is done with BIRKA in this work, thermal finite element models of a motor have been reduced by Zhou et al. [15] by capturing the most significant eigenmodes and only considering the temperature ‘hot spots’. The outline of the paper is as follows: In Section 2, BIRKA is introduced and its operation mode for large systems is explained. A norm estimate which can detect divergence in certain cases is introduced in Section 3, while two approaches for preser- ving the system’s stability are introduced in Section 4, as stability is not preserved by the algorithm. Section 5 introduces the model – an electrical motor, and the underlying thermal modelling, while Section 6 shows how the model is parameterized. Results of the reduction with BIRKA can be found in Section 7. 1.1. Notation and definitions nm kl The Kronecker product of two matrices A 2 R and B 2 R is defined as 0 1 a B ... a B 11 1m B . . C . . A#B ¼ : @ A . . aBa B n1 nm nm The vectorization operator of a matrix A 2 R is: vecðÞ A ¼½ a ... a a ... a ... a ... a 11 n1 12 n2 1m nm nn The notation λ (A),i= 1, …, n, is used for the eigenvalues of the matrix A 2 R and λ i i nn (A, B),i= 1, …, n, for the eigenvalues of the pencil (A, B) with A; B 2 R . Similarly, nm nn σ (A),i= 1, …, n, are the singular values of the matrix A 2 R . If the matrix A 2 R is positive (semi)definite, this is written as A ≥ 0or A> 0, respectively. The analogue notation holds for negative (semi)definite matrices. The spectral norm of a matrix is denoted as kk M ¼ σ ðÞ M , and the spectral norm condition number of a matrix as max σ ðÞ M 1 max κ ðÞ M ¼kk M M ¼ , where σ ðÞ M =σ ðÞ M denotes the maximal/ 2 max min 2 σ ðÞ M 2 min minimal singular value of M. 106 A. Bruns and P. Benner Definition 1.1. For a linear system (1) the H -norm is defined as: 2 T kk H ¼ tr HðÞ iω HiðÞ ω dω (6) 2π −1 where H(s)= C(sE − A) B is the transfer function corresponding to the linear system (1). Definition 1.2. For a bilinear system (5) the H -norm is defined as: Z Z 1 m 1 1 X X ðÞ l ;...l ðÞ l ;...l 2 1 k 1 k ¼ tr ... g ðg Þ ds ... ds (7) kk bil 1 k H k k 0 0 k¼1 l ;...; l ¼1 1 k As As ðÞ l ;...l 1 1 k As N e k1 N ...e k l l 1 2 with g ðÞ s ; ... ; s¼ Ce b : 1 k l k k 2. The bilinear interpolatory rational Krylov algorithm The objective of MOR is to generate a reduced model that approximates the original one with as small deviation as possible. This deviation can be measured, for example, in terms of the error in the H -norm. For a bilinear system (5), necessary H optimality conditions, 2 2 which minimize this error locally, have been derived in [11,12]. The reduced models generated with BIRKA fulfil these conditions. The reader should note that for linear, parametric systems optimality conditions are not completely known yet. A first result can be found in [4], where the condition is given using a norm on H #L ðÞ Ω , with parameter 2 2 domain Ω. The finite element discretization of industrial models leads to a system with E ≠ I . Due to their large dimension, the inversion of E would be numerically −1 expensive or even impossible and would in addition lead to full matrices E A and −1 E B. Hence, models with nonsingular E ≠ I are considered, but explicit inversion of E needs to be avoided. This can be accomplished by a modification of the original BIRKA [12] which incorporates the E in the iteration as it can be seen in Algorithm 1. There, projection matrices V and W are iteratively computed. A set of starting matrices (E = I,A ,(N ) ,B,C ) of the reduced order r is required, which is chosen randomly or r r r r k r r by an informed guess. The stopping criterion is based on the change of the eigenvalues Λ −1 of the diagonalization A = S ΛS of the iteratively generated matrices. The algorithm is stopped if the change in the eigenvalues Λ of two subsequent models is smaller than a predefined tolerance. The projection matrices for the generation of the next reduced model are calculated by T T e ~ vecðÞ V ¼ I # A  Λ# E  N # N B # B vecðÞ I ; (8a) r k k m k¼1 Mathematical and Computer Modelling of Dynamical Systems 107 T T T T e ~ vecðÞ W ¼ I # A  Λ# E  N # N C # C vec I : (8b) r k k p k¼1 Algorithm 1 [12] Bilinear IRKA for systems with E ≠ I, E nonsingular Input: E, A, N ,B,C,E,A,(N ) ,B,C k r r r k r r opt opt opt opt opt Output: E ¼ I ; A ;ðÞ N ; B ; C r r r r k r r 1: while (Change in Λ >0) do 1 T 1 1 ~ ~ ~ 2: A ¼ SΛS ; B ¼ S B ; C ¼ C SN ¼ S ðÞ N S r r r k r T T ~ ~ 3: vecðÞ V ¼ I #A  Λ#E  N #N B #B vecðÞ I r k k m k¼1 1 T T T T T 4: vecðÞ W ¼ I #A  Λ#E  N #N C #C vec I r k k p k¼1 5: V = orth(V), W = orth(W) % orth computes an orthonormal basis 1 1 1 T T T T T T 6: A ¼ðÞ W EV W AV ; ðÞ N ¼ðÞ W EV W N V ; B ¼ðÞ W EV W B; C ¼ CV r r k r r 7: end while opt opt opt opt 8: A ¼ A ;ðÞ N ¼ðÞ N ; B ¼ B ; C ¼ C r r r r r r k k r r Due to the Kronecker product, the calculation of the projection matrices V and W is not possible for large systems. In [12] an iterative method was proposed to overcome this difficulty. For the derivation a Neumann Series is used in the following way: T T ~ ~ vecðÞ V ¼ I # A  Λ# E  N # N B # B vecðÞ I r k k m k¼1 "# ! 1 m X X ¼ ðÞ I # A  Λ# E N # N r k k |{z} i¼0 k¼1 ðÞ ðÞ I # A  Λ# E B # B vecðÞ I r m ¼ðÞ I # A  Λ# E B # B vecðÞ I þ r m |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V (9) ðÞ I # A  Λ# E N # N vec V r k k k¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V j1 ...þðÞ I # A  Λ# E N # N vec V þ ... r k k k¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} vecðÞ V ¼ vec V ; j¼1 108 A. Bruns and P. Benner 1 m where (×) is only valid if ðÞ I # A  Λ# E ð N # N Þ <1 holds. In prac- r k k k¼1 tice, the infinite sum is truncated after an appropriate number of additions. The columns of the summands V are now calculated without using any Kronecker products: V ¼ðÞ λ E  A BB i i 2 1 V ¼ðÞ λ E  A N V N i k k k¼1 j 1 j1 V ¼ðÞ λ E  A N V N ; for i ¼ 1; ... ; r: i k k k¼1 This calculation can be executed in the same way for vec(W). The same projection matrices are calculated using the Truncated BIRKA proposed by Flagg [16]. The large matrices (−λ E − A) can be factorized by an LU-decomposition so that V can be efficiently calculated. For large models (n >10 ), this approach might not be feasible due to its memory demands. One can then solve the linear systems by an iterative method (cf. Saad [17]). 3. Estimation of the norm for the Kronecker product approximation Approximating the Kronecker product as in (9) can lead to divergence if ðÞ I # A  Λ# E N # N   1: r k k k¼1 It is advisable to check whether this norm remains smaller than one during the execution of BIRKA, as divergence might lead to poor reduced order models. However, a direct rnrn calculation of the norm would involve the inversion of ðI #A  Λ#EÞ2 R which might not be feasible for large systems due to memory problems. Hence, the calculation of the Kronecker product has to be avoided. To this aim, the following norm estimation is introduced: ðÞ I # A  Λ# E ð N # N Þ r k k k¼1 ðÞ I # A  Λ# E ð N # N Þ r k k k¼1 (10) ðÞ I # A  Λ# E  N # N r k k k¼1 ðÞ I # A  Λ# E N kk N : r k k |{z} k¼1 see below If the last expression is smaller than 1, the algorithm is definitely usable. We have thus derived a sufficient condition. Mathematical and Computer Modelling of Dynamical Systems 109 The normðÞ I # A  Λ# E can be calculated without explicit inversion of the matrix. The following Lemmata (cf [18]. Chapter I.4 and [1] Chapter 3) will be used to establish the new result for the calculation of the corresponding norm in Proposition 3.3. nn Lemma 3.1. For M 2 C nonsingular, M ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : min λ MM i¼1...n i Lemma 3.2. For a normal matrix M: M ¼ : min jj λðÞ M i¼1...n i By using the two Lemmata 3.1 and 3.2, we derive the following proposition, which will be used for the calculation of the norm ofðÞ I # A  Λ# E . nn Proposition 3.3. For A; E 2 R , symmetric, D = diag(d , …,d ): 1 r ðÞ I # A  D# E ¼ ; where jλ ðA  d EÞj min k for d 2 R qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi θ ¼ min : for d 2 C k¼1...r λ ððA  d EÞðA  d EÞ Þ k min k k Proof. The above matrix can be written as follows: 2 3 A  d E 0 ... 0 6 7 0 A  d E ... 0 6 7 ðÞ I #A  D#E¼ 6 . . 7; . . 4 5 . . 0 ... 0 A  d E with d 2 C.For d 2 R it is obvious that (−A − d E) is normal due to A, E symmetric, k k and thus Remark 3.2 can be used for calculating λ (−A − d E). For d 2 C the min k eigenvalue λ ðÞ A  d EðÞ A  d E is determined using Lemma 3.1. Taking the min k k minimum of all calculated eigenvalues and inverting it concludes the proof. ⊚ The calculation ofðÞ I # A  Λ# E can now be done by Proposition 3.3 using the MATLAB [19] function eigs. For the estimation of the norm as given in (10) it remains to calculate the norms of N and N , which is done in MATLAB with the functions k k normest and norm, respectively. For randomly chosen initial values the norm estimate is possibly greater than 1. However, as Λ and N change towards their optimal values, the norm estimate improves. At least two or three iterations should therefore be performed to check whether the norm for better approximations of Λ and N is smaller than 1. k 110 A. Bruns and P. Benner 4. Stability In contrast to the observations in [12], unstable systems have been encountered when applying BIRKA to industrial problems. Hence, a concept for stability preservation for the reduction with BIRKA is needed. As usual, an (asymptotically) stable linear system will be a system for which it holds Re(λ (A, E)) < 0 for i= 1, … ., n. Therefore, the parametric system in (4) is asymptotically stable if Reðλ ðA þ h N ; EÞÞ < 0 for i= 1, … ., n, and a bilinear system which was i k k k¼1 derived from a parametric system as in (5) is stable if Reðλ ðA þ u N ; EÞÞ < 0 for i k k k¼1 i= 1, … ., n. In general for bilinear systems, the inputs u can be time dependent and the definition of stability for bilinear systems will be stated in terms of the system’s input– output behaviour: Definition 4.1([11]). A bilinear system Ex_ðtÞ¼ AxðtÞþ N u ðtÞxðtÞþ BuðtÞ; k k (11) k¼1 yðtÞ¼ CxðtÞ; is said to be bounded-input-bounded-output (BIBO) stable, if for any bounded input u(t), the output y(t) is bounded on [0, ∞). For the special bilinear systems that result from parametric systems, it is possible to deduce a relation between the eigenvalues of the matrices A and A þ u N . k k k¼1 As N = 0for u resulting from the original linear inputs, only the inputs k k that are time independent will be taken into account and a comparison of the linear and bilinear eigenvalues makes sense. The next Theorem 4.2 and Corollary 4.3 are used to show Proposition 4.4 providing results for the distance between the considered eigenvalues and the stability of the bilinear system in terms of the eigenvalues. nn –1 Theorem 4.2 (Bauer–Fike [20]). If µ is an eigenvalue of M þ S 2 C and X MX = diag(λ , … ., λ ), then 1 n min jλ  μj κ ðXÞjjSjj : (12) i 2 i¼1;...;n –1 nn Corollary 4.3. Let X MX = diag(λ , … ., λ ), and M þ S 2 C . For every eigen- 1 n value λ(M + S) an eigenvalue λ (M) exists such that |λ (M) – λ(M + S)| ≤ κ (X)||S|| . i i 2 2 For a bilinear system with constant inputs u , we derive the following proposition: –1 Proposition 4.4. Let A = X diag(λ , …, λ )X with Re(λ (A)) < – c <0 for all i = 1, … ., 1 n i n. For any j 2 {1, … ., n} there exists an i 2 {1, … ., n} such that jλ ðAÞ λ ðA þ u N Þj<c i j k k k¼1 Mathematical and Computer Modelling of Dynamical Systems 111 if jjujj jjN jj < : (13) 2 2 κ ðXÞ k¼1 In addition Reðλ ðA þ u N ÞÞ <0 for j = 1, …, n. j k k k¼1 Proof. With Corollary 4.3 one concludes: m m X X jλ ðAÞ λ ðA þ u N Þj  κ ðXÞjj u N jj i j k k 2 k k k¼1 k¼1 κ ðXÞjjujj jj N jj 2 k 2 2 k¼1 <c: Assume Reðλ ðA þ u N ÞÞ ≥ 0 for one fixed j 2 {1, … ., n}. As c < |Re(λ (A))| j k k i k¼1 for all i =1, …,n and for j there exists i such that jλ ðAÞ λ ðA þ u N Þj<c one i j k k k¼1 calculates: c< jReðλ ðAÞÞj jReðλ ðAÞÞj þ Reðλ ðA þ u N ÞÞ i j k k k¼1 ¼jReðλ ðAÞÞ  Reðλ ðA þ u N ÞÞj i j k k k¼1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X X m 2 m 2 ðReðλ ðAÞ Reðλ ðA þ u N ÞÞ þðImðλ ðAÞÞ  Imðλ ðA þ u N ÞÞ i j k k i j k k k¼1 k¼1 <c; which leads to a contradiction. Therefore, Reðλ ðA þ u N ÞÞ < 0 holds. □ j k k k¼1 For systems with E= I and sufficiently small inputs u and matrices N (cf. (13)), the n k k bilinear system remains stable and every eigenvalue of the bilinear system lies in a neighbour- hood of an eigenvalue of the linear system. For E nonsingular, the Proposition 4.4 remains –1 m valid for E A and u E N . Hence, it will be assumed that the eigenvalues of k k k¼1 m –1 1 1 E A þ u E N and E A are sufficiently close and use stability preserving methods k k k¼1 for linear systems, which works quite well for the considered thermal systems. 4.1. Stability preservation using the systems Gramians For linear systems (i.e. N = 0, k= 1, …, m), stability can be preserved by using the following result due to Yousefi [21]. Basically, Villemagne and Skelton [22] already stated it, whereas Gugercin [23] used it in the context of an interpolatory approach. Yousefi incorporated the fact that the eigenvalues of the reduced model will not exceed a certain value σ. 112 A. Bruns and P. Benner Proposition 4.5 ([21]). Given a linear stable system (A, B, C) with Re(λ (A)) < – σ < 0 for nk T –1 i= 1, …, n. Then for any arbitrary full row rank matrix V 2 R and W = QV(V QV) , T T where Q = Q > 0 satisfies A Q + QA + 2σQ < 0, the reduced model (A,B,C ) is stable r r r and A =W AV satisfies Re(λ (A )) < – σ for i = 1, …, r. r i r For positive semidefinite Q, the Proposition stays valid, if one assumes V QV to be invertible. We generalize it for a system with E ≠ I, E nonsingular, positive semidefinite Q T T and Q ¼ V E QEV nonsingular, which – up to the authors’ knowledge – has not been stated elsewhere. Proposition 4.6. Given a linear stable system (E, A, B, C) with E nonsingular and Re(λ (A, E)) < – σ < 0 for i =1, …, n. Let Q = Q ≥ 0 satisfy T T T A QE þ E QA þ 2σE QE  0: (14) nk T T Then for any arbitrary full row rank matrix V 2 R with Q ¼ V E QEV nonsingular ,A,B,C ) generated with (and therefore Q>0) the reduced model (E r r r r T T 1 W ¼ QEVðV E QEVÞ is stable and satisfies Re(λ (A,E )) ≤ – σ for i =1, … ., r. i r r The proof of the Proposition follows exactly the proof of Yousefi (cf. Proposition 4.5). However, as we have introduced two generalizations – the presence of the E matrix and the nonstrict Lyapunov inequality (cf. Equation (14)) – we state it here for completeness. Proof. Multiplying Equation (14) with V and V and making use of T T 1 T T T T T T T T I ¼ðV E QEVÞ ðV E QEVÞ ¼ðV E QEVÞ ðV E QEVÞ; leads to: T T T T T T V A QEV þ V E QAV þ 2σV E QEV  0 T T 1 T T T T T T T T T T ) V A QEVðV E QEVÞ ðV E QEVÞ ðV E QEVÞ ðV E QEVÞ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ^ W T 1 T T T T T T T T T T þ V E QEVðV E QEVÞ V E QEV ðV E QEVÞ V E Q AV |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} W ^ W T T 1 T T T T T T T T T T þ 2σV E QEVðV E QEVÞ ðV E QEVÞ ðV E QEVÞ V E Q EV  0 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} W ^ W T T T T T T T T T ^ ^ ^ ) V A WQW EV þ V E WQW AV þ 2σV E WQW EV  0 T T T ^ ^ ^ ) A QE þ E QA þ 2σE QE  0 r r r r r r ^ ^ )ðA þ σE Þ QE þ E QðA þ σE Þ 0: r r r r r r Mathematical and Computer Modelling of Dynamical Systems 113 T T T T T r Using the identity E ¼ W EV ¼ðV E QEVÞ V E QEV ¼ I , let λ and v be any r r i eigenvalue and eigenvector of [A + σI ], then: r r T T ^ ^ ^ ^ ðA þ σI Þ Q þ QðA þ σI Þ 0 ) v ðA þ σI Þ Qv þ v QðA þ σI Þv  0 r r r r r r i r r i i i r r ^ ^ ) λ v Qv þ λ v Qv  0 i i i i i i r r )ðλ þ λ Þv Qv  0 i i i ) 2Reðλ Þv Qv  0 i i ðv Qv >0Þ) Reðλ Þ 0: i i Consider the eigenvalues of the reduced system (which are the eigenvalues of A as r r r E = I ) With λ v ¼ðA þ σI Þv this leads to A v ¼ λ v  σv ¼ðλ  σÞv.As r r i r r i r i i i i i i i r r Reðλ Þ 0 and – σ < 0, one can conclude that Reðλ  σÞ < 0 and therefore that the i i reduced system is stable. □ The dual result is also true: Proposition 4.7. Given a linear stable system (E, A, B, C) with E nonsingular and Re(λ (A, E)) < – σ < 0 for i = 1, …, n. Then for any arbitrary full row rank matrix nr T W 2 R and P = P ≥ 0 which satisfies T T T APE þ EPA þ 2σEPE  0 (15) T T and nonsingular P ¼ WEPE W , the reduced model (E,A,B,C ) generated with r r r r T T T T V ¼ PE W ðWEPE W Þ is stable and satisfies Re(λ (A,E )) < – σ for i = 1, …, r. i r r Proof. The proof is analogous to the one of Proposition 4.6. □ For the calculation of the projection matrix W the following Lyapunov equation is solved: T T ðA þ σEÞ QE þ E QðA þ σEÞ¼C C  0; (16) for a σ < |Re(λ (A, E))|. This leads to the calculation of the projection matrix max T T W ¼ QEVðV E QEVÞ as in Proposition 4.6. The solution of the Lyapunov Equation (16) is positive semidefinite, as the shifted system (A + σE, E) remains asymptotically stable. Equation (16) can be solved using the low rank alternating direction implicit (ADI) iteration (cf. for example [24, 25]) which generates a low rank factor Z, such that Q ≈ Z Z. The calculated low rank 114 A. Bruns and P. Benner T T T matrix Q  V E Z ZEV can be singular. This always occurs in the case where rk(Z) <rk T T 2 (V) =r. Even if rk(V) ≥ rk(Z) one cannot conclude that V E QEV is nonsingular , but for rk(V) relatively small compared to rk(Z) it is often the case. The proposed stability preservation method will be used for one of our models in Section 7.1.1. Solving large Lyapunov equations is numerically demanding. For large systems (n > 500, 000) it might be impossible – even with highly developed methods such as the ADI algorithm. Hence, this stability preserving method will reach its limitations when the systems dimensions are getting too large. 4.2. Stability preservation via one-sided projections In the special case of symmetric matrices E, A and N and positive definite E, another possibility for preserving stability is using only a single projection matrix (one-sided method). The matrices of the thermal system provided in Equation (24) in Section 6 have exactly these properties, and therefore this stability preservation is of interest. Proposition 4.8 ([26]). Given a linear system (i.e. N =0) with A, E symmetric, the system T T is asymptotically stable if and only if E = E >0 and A = A <0. nn nr Corollary 4.9 (4.4 [18],). Let A 2 R be a symmetric matrix. Let V 2 R have orthonormal columns. Let A = V AV. Then λ ðAÞ λ ðA Þ λ ðAÞ; i ¼ i; ... ; r (17) i i r iþnr T T Corollary 4.10 ([26]). Given a linear system with E = E >0 and A = A with λ (A) < 0, nr T T for i = 1, …, n. Let V 2 R have orthonormal columns. Let A =V AV and E =V EV. r r Then the reduced system is asymptotically stable. Hence for linear systems with A, E symmetric and E positive definite stability can be preserved via one-sided projections. Changing to bilinear systems, one needs the follow- ing results. nn Proposition 4.11 (Weyl [18], ). Let A; B 2 R be two symmetric matrices. Let λ (A) and λ (B) for i = 1, …, n be the eigenvalues of A and B with λ ðAÞ    λ ðAÞ and i 1 n λ ðBÞ    λ ðBÞ. Then it holds: 1 n λ ðA þ BÞ2½λ ðAÞþ λ ðBÞ; λ ðAÞþ λ ðBÞ for i ¼ 1; ... ; n (18) i i n i 1 Corollary 4.12 ([18]). Under the assumptions of Proposition 4.11 it holds jλ ðA þ BÞ λ ðAÞj  jjBjj for i ¼ 1; ... ; n (19) i i As in the context of the Bauer–Fike Theorem 4.2, the eigenvalues of a bilinear system, derived from a linear, parametric system can be related to the eigenvalues of the linear system. Using the previous results, we obtain the following Corollaries: Mathematical and Computer Modelling of Dynamical Systems 115 nn nn Corollary 4.13. Let u 2 R for k = 1, …, m, A 2 R and N 2 R symmetric with k k nr eigenvalues 0> λ ðAÞ    λ ðAÞ and λ ðN Þ λ ðN Þ. Given that V 2 R 1 n 1 k n k T T has orthonormal columns and A =V AV and (N ) =V N V, then it holds r r k k m m X X jλ ðA þ u ðN Þ Þ λ ðA Þj  jjujj jjN jj : (20) i r k r i r k k 2 2 k¼1 k¼1 Proof. Corollary 4.12 leads to m m X X jλ ðA þ u ðN Þ Þ λ ðA Þj  jj u ðN Þ jj i r k r i r k r k k 2 k¼1 k¼1 jjujj jjðN Þ jj : 2 k 2 k¼1 As N and (N ) are symmetric, they are normal and therefore fulfil ||(N ) || = max k r k r k 2 i =l, |λ ((N ) ) | = max{|λ ((N ) )|,|λ ((N ) )|} and ||N || = max |λ (N )|= max …, r i r k 1 r k r r k k 2 i =1, …,r i k {|λ (N )|, |λ (N )|}. With Corollary 4.9 one concludes λ ðN Þ λ ððN Þ Þ λ ðN Þ. This 1 k i r n k 1 k r k leads to jjðN Þ jj jjN jj and therefore Equation (20) holds. □ r k k 2 2 In addition one obtains: Corollary 4.14. Under the assumptions of Corollary 4.13 let c 2 R with c < P P m m jλ ðA Þj ¼ jλ ðA Þj. If jjujj jjN jj <c then λ ðA þ u ðN Þ Þ < 0. max r 1 r k i r k r 2 k¼1 2 k¼1 k Proof. Assume λ ðA þ u ðN Þ Þ 0 and calculate using Equation (20): i r k r k¼1 k c< jλ ðA Þj<jλ ðA Þj 1 r i r jλ ðA Þj þ λ ðA þ u ðN Þ Þ i r i r k r k¼1 ¼jλ ðA Þ λ ðA þ u ðN Þ Þj <c: i r i r k r k¼1 This leads to a contradiction and therefore λ ðA þ u ðN Þ Þ < 0. □ i r k r k¼1 k If jjujj jjN jj is sufficiently small, one can assume that the eigenvalues 2 2 k¼1 P P m m λ ðA þ u ðN Þ Þ and λ (A ) are similar. In addition, λ ðA þ u ðN Þ Þ < 0. i r k r i r i r k r k¼1 k k¼1 k Hence, starting from a model with λ (A) < 0 for all i =1, …, n and jjujj jjN jj i k 2 k¼1 2 sufficiently small leads to stable reduced bilinear systems. 4.3. Stability preservation – the workflow As the reduced models, that have been calculated with the stabilization process using the Gramians, are in most cases better than those generated by a one-sided approach, the workflow in Figure 1 applies. The reader should note that l r indicates the fact, that the matrix Q (cf. Section 4.1) can always be singular, but for the case of l r, it is more likely that Q is invertible. 116 A. Bruns and P. Benner Fix a reduced order r. Is solving equation (16) by an ADI iteration possible? Yes No Solve the equation (16) and deter- Use a one-sided approach. mine the rank rk(Z)= l. l  r l < r Stability preservation via Proposition 4.6. Figure 1. Proposed workflow for stabilization. 5. Heat transfer modelling 5.1. Mode of operation of an electrical motor An electrical motor converts electrical energy into mechanical work, which is produced by the interaction of an electrical current and a magnetic field. One part of the motor – the so- called stator – consists of several coils wound around an iron core. When a voltage is applied, a current is induced in the coil. Inside the counterpart – the so-called rotor – a magnetic field is generated either by a permanent magnet or by an electromagnet. The interaction of this magnetic field with the current in the stator results in a rotation of the rotor. See Figures 2 and 3. Actuating the motor with electrical currents leads to an increase in temperature in its different components due to thermal losses. It is important to analyse the influence of this temperature change on the materials of the motor, as it affects its lifespan. This is achieved by carrying out a thermal analysis. 5.2. Thermal modelling of an electrical motor For a thermal analysis, several physical effects have to be considered and can be modelled based on the three main types of heat transfer: heat conductance, convection and radiation. magnets coil stator rotor Figure 2. Drawing of a slice through an electrical motor. Mathematical and Computer Modelling of Dynamical Systems 117 (a) Drive unit and generator in one: the (b) Generator for commercial vehicles. The Bosch integrated motor generator. configuration of the coils is the same as for an electrical motor Figure 3. Two Bosch components illustrating the structure of an electrical motor. Pictures: Bosch. © [Bosch]. Reproduced by permission of Robert Bosch GmbH. For a broad overview of heat and mass transfer, see for example the book of Baehr and Stephan [27]. The main heat sources in the electrical motor are thermal losses, resulting from the current in the coil of the stator and/or rotor. The motor has to fulfil various operational requirements and therefore different current profiles have to be considered. The tempera- ture on certain parts of the motor (for example, the flange) should not exceed a specified upper limit, because these parts are in contact with other temperature sensitive compo- nents. This upper limit is built into the model as a fixed temperature. The motor is surrounded by air, therefore convection has to be considered. In general, convection describes the heat transfer between a fluid (here: surrounding air) and a solid. This effect can be modelled via a so-called heat transfer coefficient h. Different environ- ments have to be analysed, as the motor needs to work in a large temperature range (arctic winter, tropical summer); therefore, different ambient temperatures are examined in the model. Varying the heat transfer coefficients represent different cooling strategies or different interaction scenarios of the motor with its environment. Various parts of the motor are not directly attached and the resulting thermal resistance has to be modelled by a contact heat transfer coefficient. Varying this parameter, the small gap between the two materials can be considered as filled with air or an insulation material. The motor is built from various materials such as steel, copper and plastics. These materials have different properties, among others the density ρ, the specific heat C and the thermal conductivity k. Here, these material parameters will not be varied. As the motor temperature remains relatively small, the effect of radiation is not of great impor- tance, and therefore it will not be analysed. 6. Model parametrization To determine the thermal behaviour of the motor while in operation, the temperature field T(x, t) (dependent on location x and time t) has to be examined. The temperature field T(x, t) within a domain Ω 2 R for times t 2½ 0; t can be calculated by using the heat end equation @TxðÞ ; t ρC ¼ kTxðÞ ; t þSxðÞ ; t : (21) @t 118 A. Bruns and P. Benner using constant material properties ρ,C,k and a heat source S, which in case of an electrical motor corresponds to the thermal losses of the coil of stator or rotor. The derivation of the heat equation can be found, e.g. in [27]. On interfaces and outer surfaces, now called boundaries and denoted as Γ 2 R , different conditions have to be given, depending on the situation of interest. They are mathematically formulated as follows: ● Dirichlet boundary conditions: Tðx; tÞ¼ T ðtÞ on the boundary Γ : D D These conditions correspond to fixed temperatures on surfaces. Neumann boundary conditions: @TxðÞ ; t k ¼ q_ in Γ ; N N @n where q_ is a given heat flux on the boundary. In the considered model of the electrical motor, this kind of boundary condition does not appear. Robin boundary conditions: @TxðÞ ; t k hTðÞ  T in Γ ; 1 R @n where h denotes the heat transfer coefficient. The temperature T is the temperature of the surrounding domain. Interface conditions: A thermal resistance between two surfaces can be modelled on the interface by a thermal contact conductance coefficient h . The interface I will be considered as two surfaces: I with temperature T and I with temperature T . The following equation 1 1 2 2 applies: @TjIðÞ x; t @TjIðÞ x; t 1 2 k ¼k ¼ h TxðÞ ; tj TxðÞ ; tj : 2 1 c I I 1 2 @n @n A spatial discretization of the heat Equation (21) is then performed by using the finite element method. Incorporating the boundary and interface conditions and using a finite element basis function ψ; ψj ¼ 0, the weak formulation reads as follows: Z Z Z @TxðÞ ; t ψðÞ x ρC dx þ ψðÞ x  kTxðÞ ; t dx þ ψðÞ x hTðÞ x; t ds @t Ω Ω Γ Z Z þ ψðÞ x h TxðÞ ; t ds  ψðÞ x h TxðÞ ; t ds (22) c c I I 1 2 Z Z Z ¼ ψðÞ x SxðÞ ; t dx þ ψðÞ xðÞ qN _ ds þ ψðÞ x hT ds: Ω Γ Γ N R The material parameters ρ,C and k are considered as constant. With finite element basis functions ψ ðxÞ the temperature is approximated as follows, k Mathematical and Computer Modelling of Dynamical Systems 119 Tðx; tÞ T ðtÞψ ðxÞ: k¼1 This leads to the following discretized equation: 2 3 SðtÞ hT 6 7 ETðtÞ¼ ðA þ hN þ h N ÞTðtÞþ B: : (23) 1 c 2 4 5 q_ R R R R As ψ ðxÞψ ðxÞdx ¼ ψ ðxÞψ ðxÞdx and ψ ðxÞ ψ ðxÞdx ¼ ψ ðxÞ ψ ðxÞdx; k j j k k j j k the matrices E, A and N for the considered class of systems are symmetric and E is in addition positive definite. One should note that a thermal conductance can lead to a singular matrix A, whereas A + h N is nonsingular. Hence a new parameter h ¼ h  s is introduced, and c 2 c c ~ ~ A þ h N with A ¼ A þ sN andafixedshift s will be considered. c 2 2 7. Results The thermal analysis is carried out using COMSOL Multiphysics, version 3.5a [28]. By exporting several matrices from COMSOL and a thorough analysis of the underlying equations, it is possible to reconstruct a parametric model with variable parameters and loads. 2 3 h T 1 1 6 7 ! 6 7 q . v 6 7 X X 6 7 _ ~ ~ ETðtÞ¼ A þ h N þ ðh Þ N TðtÞþ B  ; 6 7 i i c k k h T q 1 6 7 (24) i¼1 k¼qþ1 6 7 4 5 LðtÞ yðtÞ¼ CTðtÞ; where q is the number of heat transfer coefficients h,and v – q is the number of contact heat transfer coefficients h .In Figure 4, the modelled motor part is shown. One can see parts of stator, coil, housing and some insulation parts. The following loads and parameters need to be considered. On top of the housing a temperature T is fixed to include a specified maximum temperature. The coil losses L(t) are incorporated into the coil. Heat transfer by convection is modelled at seven different locations, for example on coil and housing, resulting in 7 heat transfer coefficients (i.e. q = 7). Thermal resistance is incorporated at six different locations, for example between the insulation and the stator or the insulation and the coil (i.e. v= 13). The size of the model is n = 41, 199 and the original transient COMSOL simulation for one parameter setting takes about 90 minutes. There are two different models of the electrical motor that have been examined. The first one is simulated using only heat transfer coefficients as parameters and ignoring the effects of thermal resistance between some motor parts. This leads to a model with 7 120 A. Bruns and P. Benner Figure 4. The Comsol model simulates the heat transfer in a stator slice, without the rotor. parameters and 4 loads. The second model simulates the thermal resistances and therefore contact heat transfer coefficients have to be considered. Hence, the resulting model has 13 parameters and 4 loads. Each of the resulting parametric systems (24) is reformulated as a bilinear system by using the procedure explained in [7] and afterwards reduced using BIRKA (Algorithm 1) with the calculation of the projection matrices V and W as explained in (9). The infinite sum is truncated after 10 summands. The calculations were performed using MATLAB [19] on 12CPUs with 3GB RAM each. 7.1. Model 1 – no contact heat transfer coefficients 7.1.1. General results The first reduction has been executed using the stability preservation methods in Section 4.1. Hence, the stability of the original model is preserved by calculating the projection matrix W as proposed in Proposition 4.6. The reduction took about 11 hours. This –7 corresponds to 16 iterations and a change in the eigenvalues of less than 10 . The reduced model has order r = 50 and can be simulated within 10 seconds, which corresponds to a speed-up of over 500 compared to the original simulation time of 90 minutes. Compared to the reduction time of 11 hours, the original model could have been simulated about 8 times. As an optimal parameter set for predefined operating conditions is determined within an optimization process, it requires up to 10, 000 simulations of the model – the reduced model pays off. Now the results at four different locations will be examined. The temperatures at the front of the stator, at the coil and at two different points on the insulation have to be monitored. The original solution is compared to the solution which was computed using the reduced model. Figure 5 shows that there is only a small deviation which can be seen in the error plots. The relative error in temperature is smaller than 0.07 K, corresponding to a relative maximal error of less than 2  10 . It is important to make sure that the reduced model gives reliable results over a wide range of parameter values and inputs. The heat transfer coefficients in our model range from 5 to 100, and the coil losses L(t) and the ambient temperature T also need to be varied. For all these variations the reduced ∞ Mathematical and Computer Modelling of Dynamical Systems 121 Figure 5. Different outputs – original model compared to reduced model and relative and absolute errors. model gives an excellent approximation of the full model. Take for example Figure 6, where the behaviour of the temperature for six different heat transfer coefficients on the coil is shown. The error plots show that again the relative and absolute errors are sufficiently small. 122 A. Bruns and P. Benner Figure 6. Temperature curve for six different values (5, 25, 45, 65, 85, 100[W/m K]) of the heat transfer coefficient on the coil, and the relative and absolute errors. 7.1.2. Stability preserving – comparing the different approaches As explained in Sections 4.1 and 4.2, stability can be preserved by different procedures. Here, the following approaches will be examined: gramian BIRKA: The reduced model is calculated by using V as in Algorithm 1 with Equation (9), and the matrix W is calculated using Theorem 4.6. Mathematical and Computer Modelling of Dynamical Systems 123 BIRKA-tS: The projection matrices V and W are calculated with Algorithm 1. Stability is not preserved. Hence, in every step of the iterative process the generated reduced system is saved, and a stable system is chosen from within these systems. This stable system does not necessarily exist, and even if, it is possibly not an optimal reduced system, as it is not necessarily the final reduced system. BIRKA-oS: The reduced model is calculated with Algorithm 1. Only in the last step a one- sided projection with V is used. only V: Calculate V as in Algorithm 1. In every step of the algorithm a one-sided projection is used to calculate the reduced model. For the outputs on top of the coil and on the insulation between coil and stator, all reduced models of order 40–60 have a sufficient accuracy. Figure 7 shows the relative and absolute errors of the outputs 1 and 4 on the stator front and the insulation on the top of the stator, respectively. For the reduced order 40, the original Algorithm 1 (cf. dash-dotted line in Figure 7) performs well, whereas the other approaches are still not satisfying. However, if the order is increased, the approximation improves. With the exception of the approach of only using V within the algorithm. In Figure 8, the approaches BIRKA-oS and only V are compared for different orders. For the latter one finds, that the order has to be increased up to r = 100 to get results, that are as good as for the other methods. Figure 7. Output 1 and 4 for different reduced orders error plots. 124 A. Bruns and P. Benner Figure 8. One-sided methods. 7.2. Model 2 – contact heat transfer coefficients In the second model, thermal resistance has been taken into account. Six additional parameters (contact heat transfer coefficients h ) are incorporated into the model. They W W range from 200 up to 3600 . As mentioned in Section 6, these parameters can lead 2 2 m K m K to a singular matrix A, and a shift s needs to be introduced leading to a nonsingular matrix min max A ¼ A þ sN. For every given h 2½h ; h , one shift (the centre of the interval) is c c needed for achieving a nonsingular A. Making use of the shifts might affect the stability of Mathematical and Computer Modelling of Dynamical Systems 125 ~ ~ ~ the reduced model. For A stable, one cannot guarantee that A  sN ¼ A remains stable. r r r r The stability of A and A can be connected using the statements derived in Section 4.2 and the Bauer–Fike Theorem 4.2. For this model, the stability preservation using Theorem 4.6 will fail, because the size of a reduced model will be larger than the rank of the low rank factor in the ADI iteration. Hence, stability can only be preserved using a one-sided projection. This leads to larger reduced orders than an original BIRKA would lead to. For the reduction, the following approaches are used: ● only V: The projection matrix V is calculated as in Algorithm 1. In every step of the algorithm a one-sided projection is used to calculate the reduced model. ● BIRKA-tS: The projection matrices V and W are calculated with Algorithm 1 out of a reduced model generated by only V. Stability is not preserved. Hence, in every step of the iterative process the generated reduced system is saved, and a stable system is chosen from within these systems. This stable system does not necessarily exist and, even if, it is possibly not an optimal reduced system, as it is not necessarily the final reduced system. BIRKA-oS: The reduced model is calculated with Algorithm 1 out of a reduced model generated by only V. Only in the last step a one-sided projection with V is used. The reduction was performed using the one-sided approach only V and took about 3 days and 3 hours. The reduced model has order r = 600 and can be simulated within 60 seconds, which corresponds to a speed-up of about 90 compared to the original simulation time of 90 minutes. This reduced model leads to a good approximation of the original model within the whole parameter range. This can be seen for instance in Figure 9, where the variation of the heat transfer coefficient on the coil is shown. The two approaches BIRKA-tS and BIRKA-oS use the reduced model calculated with only V and reduce it again. This leads to models of order r = 300, and additional 12 hours for the reduction process. These models can be simulated in 15 seconds, which corresponds to a speed-up of over 300 compared to the original simulation time. Figure 10 shows the original and the reduced models generated by using the above approaches, and the largest error which occurs in output 3. The reduced models generated with only V and BIRKA- tS show sufficient accuracy, whereas the approximation with BIRKA-oS is still not accurate enough. 8. Discussion and conclusion While applying BIRKA to thermal industrial models, several problems were encoun- tered. First, approximating the original calculation of the projection matrices as shown in (9) could lead to a divergence of the algorithm. This divergence can be detected by the procedure shown in Section 3; however, there still might be other causes for divergence. Second, one class of physical parameters (contact heat transfer coefficients) leads to a singular stiffness matrix A. Hence, the reduction is performed using a shifted matrix A, as explained in Section 6. Third, stability is not preserved by the algorithm and two main strategies to stabilize the reduced models are presented in Sections 4.1 and 4.2. Considering the speed-up and the accuracy of the reduced model, BIRKA seems to be an option for the reduction of industrial thermal models. As the reduction times are relatively long (especially for the second model, cf. Section 7.2), the reduction is only reasonable if a large amount of test cases with different parameter samples need to be simulated. Other methods for the reduction of parametric models use a different approach. In general, models are reduced in several parameter points, and an 126 A. Bruns and P. Benner Figure 9. Temperature curve for six different values (5, 25, 45, 65, 85, 100[W/m2K]) of the heat transfer coefficient on the coil, and the relative and absolute errors. interpolation procedure is used to obtain reduced models in new parameter points (as for example in [2–6, 29]). For models with a small parameter range, these methods might lead to smaller reduced orders, as BIRKA offers a good accuracy in a large parameter space, which might not be needed for the considered model. In addition, BIRKA is restricted to the special parameter dependence given in Equation (4). Linear models with a parameter dependence in the matrices E, B or C or with a nonaffine parameter dependence cannot be reduced directly with this algorithm, and hence other methods need to be applied. The methods given in [2–6, 29]are basedon linear Mathematical and Computer Modelling of Dynamical Systems 127 Figure 10. Reduction of model 2 with different approaches, largest error – output 3. reduction in certain parameter points. If the reduction in one sample point is not time- consuming and a small number of sampling points is needed, it is likely that the method outperforms BIRKA in terms of reduction time. However, these methods can only guarantee a small reduction error in the neighbourhood of the sampling points. BIRKA, in contrast, will yield a model with a globally small reduction error. 128 A. Bruns and P. Benner Notes 1.kk M # M ¼kk M kk M ½ 30 , Corollary 13.11. 1 2 1 2 2 2 2 0 1 0 1 0 1 1 100 00 B C T T T @ AB C @ A 2. Take n ¼ 4; E ¼ I ; V E Z ¼ 0 011 ¼ 00 @ A 0 1 1 000 10 References [1] A.C. Antoulas, Approximation of Large-Scale Dynamical Systems, SIAM, Philadelphia, 2005. [2] D.S. Weile, E. Michielssen, E. Grimme, and K. Gallivan, A method for generating rational interpolant reduced order models of two-parameter linear systems, Appl. Math. Letters. 12 (1999), pp. 93–102. doi:10.1016/S0893-9659(99)00063-4 [3] L. Daniel, O. Siong, L.S. Chay, K.H. Lee, J. White, and J. White, A multiparameter moment matching model reduction approach for generating geometrically parameterized interconnect performance models, IEEE Trans. Comput.-Aided Design Integr. Circuits Syst. 23 (2004), pp. 678–693. doi:10.1109/TCAD.2004.826583 [4] U. Baur, C.A. Beattie, P. Benner, and S. Gugercin, Interpolatory projection methods for parameterized model reduction, SIAM J. Sci. Comp. 33 (2011), pp. 2489–2518. doi:10.1137/090776925 [5] H. Panzer, J. Mohring, R. Eid, and B. Lohmann, Parametric model order reduction by matrix interpolation,at – Automatisierungstechnik. 58 (2010), pp. 475–484. doi:10.1524/ auto.2010.0863 [6] D. Amsallem and C. Farhat, An online method for interpolating linear parametric reduced- order models, SIAM J. Sci. Comp. 33 (2011), pp. 2169–2198. doi:10.1137/100813051 [7] P. Benner and T. Breiten, On H -model reduction of linear parameter-varying systems,in Proceedings of Applied Mathematics and Mechanics, WILEY-VCH Verlag, Weinheim, 2011, pp. 805–806. doi:http://dx.doi.org/10.1002/pamm.201110391 [8] S. Al-Baiyat and M. Bettayeb, A New Model Reduction Scheme for k-Power Bilinear Systems, Proceedings of the 32nd IEEE Conference on Decision and Control Vol. 1, 1993, pp. 22–27. [9] Z. Bai and D. Skoogh, A projection method for model reduction of bilinear dynamical systems, Linear Algebra Appl. 415 (2006), pp. 406–425. doi:10.1016/j.laa.2005.04.032 [10] T. Breiten and T. Damm, Krylov subspace methods for model order reduction of bilinear control systems, Syst. Control Lett. 59 (2010), pp. 443–450. doi:10.1016/j. sysconle.2010.06.003 [11] L. Zhang and J. Lam, On H model reduction of bilinear systems, Automatica. 38 (2002), pp. 205–216. doi:10.1016/S0005-1098(01)00204-7 [12] P. Benner and T. Breiten, Interpolation-based H -model reduction of bilinear control systems, SIAM J. Matrix Anal. Appl. 33 (2012), pp. 859–885. doi:10.1137/110836742 [13] Y. Bertin, E. Videcoq, S. Thieblin, and D. Petit, Thermal behavior of an electrical motor through a reduced model, IEEE T. Energy Conver. 15 (2000), pp. 129–134. doi:10.1109/ 60.866989 [14] Z. Gao, R. Colby, T. Habetler, and R. Harley, A model reduction perspective on thermal models for induction machine overload relays, IEEE Trans. Ind. Electronics. 55 (2008), pp. 3525–3534. doi:10.1109/TIE.2008.926772 [15] K. Zhou, J. Pries, H. Hofmann, Y. Kim, T.-K. Lee, and Z. Filipi, Computationally-efficient finite-element-based thermal models of electric machines, Vehicle Power and Propulsion Conference (VPPC), 2011 IEEE, Chicago, IL, IEEE, 6–9 September 2011, 1–6. Available at http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6043205&isnumber=6042961 (Accessed 2 June 2014). [16] G.M. Flagg, Interpolation Methods for the Model Reduction of Bilinear Systems, Ph.D. thesis, Virginia Polytechnic Institute and State University, Blacksburg, VI, 2012. [17] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, 2003. [18] G.W. Stewart and J.G. Sun, Matrix Perturbation Theory (Computer Science and Scientific Computing), Academic Press, San Diego, CA, 1990. [19] MATLAB, version 7.10.0 (R2010a), The MathWorks Inc., Natick, Massachusetts, 2010, ©The MathWorks, Inc. MATLAB and Simulink are registered trademarks of The MathWorks, Inc. See www.mathworks.com/trademarks for a list of additional trademarks. Mathematical and Computer Modelling of Dynamical Systems 129 Other product or brand names may be trademarks or registered trademarks of their respective holders. [20] G. Golub and C.V. Loan, Matrix Computations, 3rd ed. John Hopkins University Press, Baltimore, MD, 1996. [21] A. Yousefi, Preserving stability in model and controller reduction with application to embedded systems, Ph.D. thesis, Technische Universität München, München, 2006. [22] C.D. Villemagne and R.E. Skelton, Model reductions using a projection formulation,J. Control. 46 (1987), pp. 2141–2169. doi:10.1080/00207178708934040. [23] S. Gugercin, An iterative SVD-Krylov based method for model reduction of large-scale dynamical systems, Linear Algebra Appl. 428 (2008), pp. 1964–1986. doi:10.1016/j.laa.2007.10.041. [24] P. Benner, J.R. Li, and T. Penzl, Numerical solution of large-scale lyapunov equations, riccati equations, and linear-quadratic optimal control problems, Numer. Linear Algebra Appl. 15 (2008), pp. 755–777. doi:10.1002/nla.622. [25] J. Saak, Efficient numerical solution of large scale algebraic matrix equations in PDE control and model order reduction, Ph.D. thesis, Technische Universität, Chemnitz, 2009. [26] R. Castan Selga, The matrix measure framework for projection-based model order reduction, Ph.D. thesis, Technische Universität, München, 2011. [27] H.D. Baehr and K. Stephan, Heat and Mass Transfer, Springer, Berlin, 2006. [28] COMSOL-Multiphysics, version 3.5a, COMSOL is a registered trademark of COMSOL AB., [29] P. Benner and L. Feng, A robust algorithm for parametric model order reduction based on implicit moment matching, Reduced Order Methods for Modeling and Computational Reduction, in MS&A, A. Quarteroni and G. Rozza, eds., Springer, 2014, vol. 9, 159–186. [30] A.J. Laub, Matrix Analysis for Scientists and Engineers, SIAM, Philadelphia, 2005.

Journal

Mathematical and Computer Modelling of Dynamical Systems – Taylor & Francis

Published: Mar 4, 2015

Keywords: parametric model order reduction; stability preservation; thermal heat transfer; finite element modelling; 34K17; 93A15

References

  • A method for generating rational interpolant reduced order models of two-parameter linear systems
  • A multiparameter moment matching model reduction approach for generating geometrically parameterized interconnect performance models
  • Interpolatory projection methods for parameterized model reduction
  • Parametric model order reduction by matrix interpolation
  • An online method for interpolating linear parametric reduced-order models
  • A projection method for model reduction of bilinear dynamical systems
  • Krylov subspace methods for model order reduction of bilinear control systems
  • On H2 model reduction of bilinear systems
  • Interpolation-based H2-model reduction of bilinear control systems
  • Thermal behavior of an electrical motor through a reduced model
  • A model reduction perspective on thermal models for induction machine overload relays
  • Model reductions using a projection formulation
  • An iterative SVD-Krylov based method for model reduction of large-scale dynamical systems
  • Numerical solution of large-scale lyapunov equations, riccati equations, and linear-quadratic optimal control problems