A Block Arnoldi Algorithm Based Reduced-Order Model Applied to Large-Scale Algebraic Equations of a 3-D Field Problem

Featured Application: This work has developed an adaptive multipoint model reduction model based on the Arnoldi algorithm to obtain reduced-order models of a 3-D temperature ﬁeld. Abstract: In simulations of three-dimensional transient physics ﬁlled through a numerical approach, the order of the equation set of high-ﬁdelity models is extremely high. To eliminate the large dimension of equations, a model order reduction (MOR) technique is introduced. In the existing MOR methods, the block Arnoldi algorithm-based MOR method is numerically stable, achieving a passively reduced order model. Nevertheless, this method performs poorly when it is applied to very wide-frequency transients. To eliminate this deﬁciency, multipoint MOR methods are emerging. However, it is hard to directly apply an existing multipoint MOR method to a 3-D transient ﬁeld equation set. The implementation issues in a reduction process (such as the selection of expansion points, the number of moments matched at a point and the error bound) have not been explored in detail. In this respect, an adaptive multipoint model reduction model based on the Arnoldi algorithm is proposed to obtain the reduced-order models of a 3-D temperature ﬁeld. The originality of this study is the proposal of a novel adaptive algorithm for selecting expansion points, matching moments automatically, using a posterior-error estimator based on temperature response coupled with a network topological method (NTM). The computational efﬁciency and accuracy of the proposed method are evaluated by the numerical results from solving the temperature ﬁeld of a prototype insulated-gate bipolar transistor (IGBT).


Introduction
In computations of three-dimensional transient fields of a high-fidelity model, the order of the equation set is extremely high. It is thus infeasible to directly use the high-fidelity model to simulate the transient behavior of an electromagnetic device over a long-time span. To compromise the computational time and model accuracy performances, a model order reduction (MOR) technique is used to generate a reduced-order model to substitute the original high-fidelity one. The reduced-order model diminishes the computational complexity, utilizing less computational resources compared to the original system. Different MOR techniques have been proposed and applied to combat the complexity of interconnected problems [1][2][3] and finite-element equations of electromagnetic problems [4][5][6][7]. However, there are few works oriented for 3-D temperature fields. Moreover, in addition to the solution efficiency, an MOR approach should also be compromised of different implementation criteria. Implementation issues include the number and generating mechanism of the expansion points, the approximate order of the model for each expansion point, and so on. In this direction, some MOR techniques have already been reported. For example, the Asymptotic Waveform Evaluation (AWE), as reported in [8], has focused on solving the number and generating mechanism of the expansion points by considering complex frequency hopping (CFH). Authors in [9][10][11] have focused on the error bound and selection of expansion points for Páde-via-Lanczos (PVL). In addition, authors in [12] have focused on the error bound and selection of expansion points for Efficient Nodal Order Reduction (ENOR) by embedding proper orthogonal decomposition (POD).
However, both AWE and ENOR use a recursive expression to explicitly compute the moments of the high-fidelity model. Therefore, AWE and ENOR are not numerically stable [13]. Although PVL is a powerful and efficient algorithm, it cannot preserve some physically feasible conditions such as model passivity [13]. On the other hand, in the MOR methods based on moment matching, the Passive Reduced-order Interconnect Macromodeling Algorithms (PRIMA), based on the Arnoldi method, are numerically stable and can preserve the passivity of the high-fidelity model in RLC circuits. Nevertheless, this method is difficult to directly employ on a high-fidelity model of a three-dimensional transient field since this methodology is intended to deal with interconnect circuit equations, which are different in form from 3-D transient field equations. Furthermore, the difficulty of applications in broadband situations for PRIMA has not been explored in detail in the literature. In this regard, we propose an adaptive multipoint model reduction model based on the Arnoldi algorithm to obtain the reduced-order models of a 3-D temperature field.

Three-Dimensional Transient Temperature Field Problems
In general, the general equation set of a transient temperature filed can be formulated as [12]: where, B ∈ n×n , T ∈ n×1 , A ∈ n×n , P ∈ n×N , J ∈ N×1 , T| t=0 = 0, and their definitions are referred to [12] for space limitations. Applying Laplace transform to both sides of (1), one has (Bs + A)T(s) = PJ(s)

Moment Matching Method for 3-D Transient Field Problems
Moment matching is a projection-based method. In a moment-matching method, a subspace basis V ∈ n×r is first determined to approximate the residing manifold of the state vector, with the Petrov-Galerkin projection then used to construct the ROM. To preserve the passivity of the high-fidelity model, W = V [14] is selected, yielding a reduced-order model: where, B = V T BV, A = V T AV, P = V T P, T = TV T and r ≤ n is the size of the reduced model (3). The matrix V is computed from the following procedures. The transfer function (matrix) is One expands H(s) at point s 0 to where, M i = (−(Bs 0 + A) −1 A) i (Bs 0 + A) −1 P, i (0, 1, 2, . . .) is the ith moment of the system at the expansion point s 0 . The subspace spanned by the columns of V is range{V} = span P (s 0 ), ((Â(s 0 )) −1 A)P(s 0 ), . . . , ((Â(s 0 )) −1 A) qP (s 0 ) where,P(s 0 ) = (Bs 0 + A) −1 P,Â(s 0 ) = Bs 0 + A. In this paper, ((Â(s 0 )) −1 A) iP (s 0 ) is the ith order moment vector. From [14], one has where, H r (s) is the transfer function of the reduced model. It implies that a s 0 will guarantee that the reduced-order model achieves a relatively high accuracy of a specified frequency, with a larger q yielding a more accurate H r (s). From (7), it is obvious that H r (s) depends also on s 0 . For a fixed q, the error of H r at s is inversely proportional to the distance of s and s 0 . In other words, if s 0 is put at zero, then H r (s) is probably not accurate at high frequencies, which is explained in detail in Figure 1. where, In this paper, From [14], one has where, () r s H is the transfer function of the reduced model. It implies that a s0 will guarantee that the reduced-order model achieves a relatively high accuracy of a specified frequency, with a larger q yielding a more accurate

Multipoint Expansion
The development of a feasible methodology for a wide-frequency band ROM requires compromising the accuracy of the constructed projection basis and computational cost of the procedure.
Nevertheless, a basic MOR algorithm generally employs one expansion point in developing the projection basis. However, as aforementioned, the single-point momentmatching reduced model is a locally optimal one, and it is applicable to simulating the transient response only around the frequency in question. As the frequency band of interest increases, a large number of moments need to be included to guarantee the model's accuracy in the frequency band in question, increasing the computational and memory costs. Moreover, such a reduced model is inaccurate for simulating the transient performance of a frequency whose distance to the frequency expansion point is large. To address the locally supported approximation characteristics of a single-point moment-matching

Multipoint Expansion
The development of a feasible methodology for a wide-frequency band ROM requires compromising the accuracy of the constructed projection basis and computational cost of the procedure.
Nevertheless, a basic MOR algorithm generally employs one expansion point in developing the projection basis. However, as aforementioned, the single-point momentmatching reduced model is a locally optimal one, and it is applicable to simulating the transient response only around the frequency in question. As the frequency band of interest increases, a large number of moments need to be included to guarantee the model's accuracy in the frequency band in question, increasing the computational and memory costs. Moreover, such a reduced model is inaccurate for simulating the transient performance of a frequency whose distance to the frequency expansion point is large. To address the locally supported approximation characteristics of a single-point moment-matching approach, the Reduced-Basis Method (RBM) [15] or Balanced Truncation Proper Orthogonal Decomposition (BT-POD) [16], based on solution snapshots [17,18] from different frequency points, is introduced. However, it involves computing snapshots on a large number of points. Moreover, for each point, a factorization operation on the system matrix is required. Consequently, such an approach considerably increases computational cost.
To take advantage of both a snapshot-based and a single-point moment-matching approach, a multipoint MOR method is developed [19][20][21][22]. A multipoint MOR method obtains the projection basis from the moments of the transfer function at multiple frequency points, as explained in Figure 2.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 4 of 16 approach, the Reduced-Basis Method (RBM) [15] or Balanced Truncation Proper Orthogonal Decomposition (BT-POD) [16], based on solution snapshots [17,18] from different frequency points, is introduced. However, it involves computing snapshots on a large number of points. Moreover, for each point, a factorization operation on the system matrix is required. Consequently, such an approach considerably increases computational cost. To take advantage of both a snapshot-based and a single-point moment-matching approach, a multipoint MOR method is developed [19][20][21][22]. A multipoint MOR method obtains the projection basis from the moments of the transfer function at multiple frequency points, as explained in Figure 2.
accurate for a wide frequency range In [19], a fully adaptive scheme based on Krylov space was presented to reduce the size of a large-scale integrated circuit. This approach provides an interesting way of adaptively choosing expansion points by a posterior error estimation based on admittance or impedance. However, the computational process is rather complex and inefficient. Moreover, in 3-D transient field computations, the concern is more the precision or accuracy of In [19], a fully adaptive scheme based on Krylov space was presented to reduce the size of a large-scale integrated circuit. This approach provides an interesting way of adaptively choosing expansion points by a posterior error estimation based on admittance or impedance. However, the computational process is rather complex and inefficient. Moreover, in 3-D transient field computations, the concern is more the precision or accuracy of the solution rather than the admittance or impedance parameters. Consequently, a posterior-error estimate, based on admittance or impedance, is not suitable. In [20], a greedy multipoint MOR technique chose expansion points from an arithmetic midpoint, using a posterior-error estimation based on scattering parameters, to reduce the size of a second-order electromagnetic system. In [21], an adaptive multipoint-MOR scheme, using eigen values as the stopping criterion and an error estimator, chose the new expansion point by computing error estimators at three points every time to reduce the computational cost of large-scale inductive partial element equivalent circuits (PEEC). In [22], a multipoint MOR scheme was used to reduce the size of the finite-element equation of second-order systems. The location of the expansion point and size of the projection basis were determined based on an error estimator. However, the scheme cannot be applied directly to the transient temperature problem.
In this regard, we propose a novel block Arnoldi-based multipoint-MOR scheme for 3-D transient field problems. The proposed approach not only benefits from the work of [20][21][22], but also explores the application of Krylov subspaces. The originality of this study is a novel adaptive algorithm for determining expansion points, matching moments automatically, as well as a posterior-error estimator based on the temperature response, coupled with a network topological method (NTM).
It should be noted that the optimal expansion point should be the one that gives the largest error between the current ROM and the high-fidelity model. However, finding such an ideal point is neither practical nor feasible [9]. On the other hand, although the performance of the ROM is related to the expansion point, the Krylov-based MOR is relatively insensitive to it. It is usually chosen at the point in the right-half plane whose distance to the imaginary axis is of the same magnitude as the frequency of interest [13]. From this point of view, in our application of 3-D transient field problems, the new expansion point is determined using a logarithmic scale to minimize large errors overall.

Adaptive Scheme
To start, because of the aforementioned reasons explained in 2.3, to choose the new expansion point, the logarithm midpoint is obtained from the two old expansion points as where, f s mid denotes the midpoint frequency in logarithmic scale, with f s min and f s max denoting the minimum and maximum frequencies of the two old points, respectively. Second, the "local error" is defined as the error between the reduced-order model (ROM) and the full model at frequency, with the "global error" as the error between the ROM and the full model over the frequency interval [ f min , f max ]. The tolerances of the "local error" and "global errors" are initially set as tol l and tol g , respectively. Moreover, one defines a posteriori-error estimator, as given in (9), to determine whether the current reduced model satisfies the accuracy requirements over a certain frequency band or only at one frequency point. Moreover, where, T j (t) and T r j (t) denote the time-domain responses of the jth variable of full-and reduced-order models under an excitation signal, respectively. Parameter N moments max is defined as the limit of the number of moments calculated in one expansion point, and parameter N s max as the limit of the number of total expansion points used in the MOR procedure. The larger the number of moments at one frequency, the smaller the error between the ROM and the full model at this frequency is. Therefore, once the error of a frequency is smaller than the predefined tolerance, the algorithm stops moment matching at this frequency. In other words, N moments max is equal to the smallest number of moments under which the error between the ROM and the full model is smaller than the tolerance at the frequency. Moreover, if the actual number of the expansion points used is larger than N s max , the tolerance of the "local error" (tol l ) is reduced. To facilitate the implementation of the proposed scheme, it is explained step by step as follows.
An expansion point is first chosen (here we choose s 1 = f min ). The moment-matching procedure at s 1 is then activated sequentially until the local error between the ROM and the full model at frequency f min is smaller than the predefined tolerance. The global error (E over the frequency interval [ f min , f max ]) checking procedure is then followed. If the error is E < tol g , the algorithm stops and outputs the reduced model; otherwise f max is selected as the second expansion point. The moment-matching producer is repeated for f max until the "local error" between the ROM and the high fidelity model is smaller than the tolerance, and the subsequent block moments are aggregated to the basis formed by the previous ones, with an orthogonal normalization then performed to compress the projection basis. The global error E, between the ROM and high-fidelity model over the frequency interval [ f min , f max ], is then checked. If the error is E < tol, the algorithm stops and outputs the reduced-model; otherwise a middle frequency point f mid is determined by (8), between f min and f max selected as the third expansion point. For convenience of explanation, the serial numbers of 1 , 3 , 2 are used to represent f min , f mid , and f max , respectively. Since the moments at f min , f mid , and f max have already been added to the subspace of the projection basis, the next expansion point added to the subspace is chosen at the middle point between 1 , 3 ( f min , f mid ) by Equation (8), with the serial number of the new point marked as 1 -3 . If the transformation matrix constructed by 1 , 1 -3 , 3 , and 2 does not satisfy the accuracy tolerance, a new point is chosen, i.e., the middle point between 3 and 2 , with the serial number of the new point marked as 3 -2 . It can be seen that the middle point between every two adjacent points is used, with the serial numbers corresponding to the points as { 1 , 1 -3 , 3 , 3 -2 , 2 }. If the error is still larger than tol g , the middle point between every two adjacent points is chosen, i.e., the middle point between 1 and 1 -3 , 1 -3 and 3 , 3 and 3 -2 , and 3 -2 and 2 . The iterative procedure is repeated until the error is smaller than the predefined value.
Moreover, an orthogonalization operation is applied to the so-far generated projection basis to minimize its size by discarding its redundant vectors. Although the orthogonalization and compression procedures increase computational cost, a significant improvement in performances is obtained in the projection size and reliability of the whole reduction process. where, s 1 , · · · , s k are the expansion points used to obtain the ROM, and q 1 , · · · q k are the number of moments matched at the corresponding expansion point.
Get T from Equation (1)

Numerical Application
The proposed scheme and model are used to compute the transient temperature fields of a prototype insulated-gate bipolar transistor (IGBT) module (Figure 3a) to validate its advantages. The temperature field is computed by the network topological approach [23]. Figure 3b shows meshes of the NTM for the prototype IGBT module. In the numerical implementation of NTM, the numbers of the total elements and nodes are, respectively, 3364 and 4303, with the number of the NTM equations equal to 4303. Obviously, the solution of such a large equation set will consume an enormous amount of computational resources.  To simulate the wide-frequency-band transient response, the frequency broadband of excitation signals consisted of 100 liner-spaced samples from 10 Hz to 100 kHz. For performance comparisons, the NTM, the proposed reduced model and the conventional block-Arnoldi reduced model are used to compute the transient temperature field of this case study under the predefined, tested broadband excitation signal. In order to take both the minimum and maximum frequencies into consideration, the time step should be small enough and the total excitation time large enough. Consequently, the total excitation time was set as 0.2 s, the discrete number as 80,000 and the time step as 2.5 us. The equation was solved by ode15 s of MATLAB. In the solver settings, the relative tolerance (RelTol) was set as 10 −3 and the absolute tolerance (AbsTol) as 10 −6 . The final expansion points used were: 10, 10 5 , 10 3 , and 10 2 . As explained in the adaptive scheme, moments max N is equal to the smallest number of moments under which the error between the ROM and the full model is smaller than the tolerance at the frequency. moments max N , at 10, 10 5 , 10 3 , and 10 2 , is finally set as 6, 6, 7, and 9, respectively. Finally, max s N is set as 100 to limit the number of total expansion points used in the MOR procedure. Moreover, the local error tolerance tol l is set as 0.005, with the global error tolerance set as 0.01. Typical observing points (points 1, 2, and 3, see Figure 3b) are used for performance comparison of different ROMs. The numerical results are presented in Figures 4-9.  Tables 1 and 2. Table 1 compares the transient response of different observing points using the proposed and NTM models. Table 2 compares the absolute errors on different observing points using the proposed and conventional schemes. The numerical results in Tables 1 and 2 also validate the conclusion obtained from Figures 4-9. To simulate the wide-frequency-band transient response, the frequency broadband of excitation signals consisted of 100 liner-spaced samples from 10 Hz to 100 kHz. For performance comparisons, the NTM, the proposed reduced model and the conventional block-Arnoldi reduced model are used to compute the transient temperature field of this case study under the predefined, tested broadband excitation signal. In order to take both the minimum and maximum frequencies into consideration, the time step should be small enough and the total excitation time large enough. Consequently, the total excitation time was set as 0.2 s, the discrete number as 80,000 and the time step as 2.5 us. The equation was solved by ode15 s of MATLAB. In the solver settings, the relative tolerance (RelTol) was set as 10 −3 and the absolute tolerance (AbsTol) as 10 −6 . The final expansion points used were: 10, 10 5 , 10 3 , and 10 2 . As explained in the adaptive scheme, N moments max is equal to the smallest number of moments under which the error between the ROM and the full model is smaller than the tolerance at the frequency. N moments max , at 10, 10 5 , 10 3 , and 10 2 , is finally set as 6, 6, 7, and 9, respectively. Finally, N s max is set as 100 to limit the number of total expansion points used in the MOR procedure. Moreover, the local error tolerance tol l is set as 0.005, with the global error tolerance set as 0.01. Typical observing points (points 1, 2, and 3, see  Tables 1 and 2. Table 1 compares the transient response of different observing points using the proposed and NTM models. Table 2 compares the absolute errors on different observing points using the proposed and conventional schemes. The numerical results in Tables 1 and 2 also validate the conclusion obtained from Figures 4-9.  Table 3 gives the relative errors (norm-2 error) computed by Equation (9) on points 1, 2, and 3 of the results of the NTM model, those of the proposed reduced model (size = 28), those of the conventional block-Arnoldi reduced model (size = 20) and those of the conventional block-Arnoldi reduced model (size = 30). From the numerical results of the two conventional ROMs, it is easy to observe that a high-order reduced model produces a more accurate result compared to a low-order one. However, the high-order conventional reduced model is less accurate than the proposed one, which validates the merits of the proposed ROM scheme.
Moreover  Table 3 gives the relative errors (norm-2 error) computed by Equation (9) on points 1, 2, and 3 of the results of the NTM model, those of the proposed reduced model (size = 28), those of the conventional block-Arnoldi reduced model (size = 20) and those of the conventional block-Arnoldi reduced model (size = 30). From the numerical results of the two conventional ROMs, it is easy to observe that a high-order reduced model produces a more accurate result compared to a low-order one. However, the high-order conventional reduced model is less accurate than the proposed one, which validates the merits of the proposed ROM scheme. Moreover

Conclusions
To increase the computational speed of a high-fidelity model in solving a complex 3-D transient field problem by a numerical approach, a block Arnoldi-based multipoint-model order-reduction scheme is presented. The proposed methodology was used to reduce the size of an NTM high-fidelity model in solving the three-dimensional temperature field of a case study. The computational results have demonstrated: (1) The model reduced from the proposed methodology could achieve a higher accuracy compared to one reduced from conventional block Arnoldi in a wideband frequency. It should be noted that an MOR method is generally applicable to a linear system. However, for a nonlinear system, a linearization technique must first be implemented in the system before the application of an MOR. Additionally, once the proposed ROM has been learned in the offline stage, it can be used to simulate a transient response of the system for any excitation J(t), as long as its frequency spectrum stays in the learned-frequency interval.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.