Reduced-order modelling of parametric systems via interpolation of heterogeneous surrogates

This paper studies parametric reduced-order modeling via the interpolation of linear multiple-input multiple-output reduced-order, or, more general, surrogate models in the frequency domain. It shows that realization plays a central role and two methods based on different realizations are proposed. Interpolation of reduced-order models in the Loewner representation is equivalent to interpolating the corresponding frequency response functions. Interpolation of reduced-order models in the (real) pole-residue representation is equivalent to interpolating the positions and residues of the poles. The latter pole-matching approach proves to be more natural in approximating system dynamics. Numerical results demonstrate the high efficiency and wide applicability of the pole-matching method. It is shown to be efficient in interpolating surrogate models built by several different methods, including the balanced truncation method, the Krylov method, the Loewner framework, and a system identification method. It is even capable of interpolating a reduced-order model built by a projection-based method and a reduced-order model built by a data-driven method. Its other merits include low computational cost, small size of the parametric reduced-order model, relative insensitivity to the dimension of the parameter space, and capability of dealing with complicated parameter dependence.


Introduction
Model order reduction (MOR), as a flourishing computational technique to accelerate simulation-based system analysis in the last decades, has been applied successfully to high-dimensional models arising from various fields such as circuit simulations [1][2][3], (vibro) acoustics [4], design of microelectromechanical systems [5], and chromatography in chemical engineering [6]. The objective of MOR is to compute a reduced-order model (ROM) of small size k that captures some important characteristics of the original highdimensional model of size n (normally k n), such as dominant moments and leading Hankel singular values. When parametric studies are performed, it is desired to build a parametric ROM (pROM), which not only incorporates the modeling parameters as free parameters, but also approximates the input-output behavior of the full-order model (FOM) well for any parameter value within the domain of interest. Therefore, many parametric MOR (PMOR) methods have been proposed; see, e.g., the recent survey [7].
In general, PMOR methods can be categorized into two general types according to PMOR review paper [7]: 1. Building a single pROM by using a global basis over the parameter space. These methods [8] have received intensive research attention and proved to be efficient in many applications in the past few years, most notably the reduced basis method [9]. Despite its success in many application fields, these methods normally meet difficulty in dealing with many parameters since both the computational cost and the size of the basis matrices may grow exponentially with the dimension of the parameter space because of the curse of dimensionality. 2. Building an interpolatory pROM by interpolating local matrices or bases at parameter samples, e.g., each of which is obtained by applying nonparametric MOR at the corresponding parameters. These methods are less studied and their performance is not so satisfactory. According to [7], this approach can be again categorized into three types: a. Interpolation among local basis matrices [10]. Since the basis matrices V i ∈ R n×k are elements on the Stiefel manifold, this method interpolates V i along the geodesic on the Stiefel manifold with the help of the tangent space. However, assuming that the basis matrix evolves on the geodesic is only heuristic: it is only one of the infinite possible evolution paths on the manifold. According to our experience in numerical tests, the resulting ROM often diverges. Another disadvantage is that, this method requires the storage of all the basis matrices, which may not only be expensive, but also infeasible in many cases. First, when the ROMs are built by a non-projection-based MOR method, e.g., data-driven MOR methods or MOR methods based on physical reasoning, no bases are computed. In addition, in many practical applications, the parametric FOM is not available and what we can obtain are nonparametric FOMs built for some parameter samples, which may cause difficulty for these methods. For example, these FOMs may be obtained from discretization of partial differential equations with different meshes, may also be of different dimensions, and their realizations may be not consistent. b. Interpolation among local reduced-order state-space matrices [11,12]. As is shown in [11], directly interpolating the local reduced-order state-space matrices does not work in general. A major difficulty is that, a dynamical system has infinitely many realizations and interpolating between different realizations can result in completely wrong results. For example, a system with the same inputoutput behavior can be obtained by rearranging the rows of the matrices, but interpolating two such system matrices normally makes no physical sense, e.g., interpolating between K (p 1 ) C 0 I and 0 I K (p 2 ) C , where K (p), C are square matrices with the same size and I and 0 are corresponding identity and zero matrices, respectively. Therefore, a natural idea is to apply a congruence transformation to obtain consistent bases, i.e., by solving an optimization problem to get the transformation matrices, and then interpolate these consistent ROMs on some manifold [12]. Another choice is to conduct a singular value decomposition (SVD) on the union of all basis matrices to calculate the dominant "global subspace", onto which we re-project all ROMs and conduct interpolation [11]. However, like the methods discussed in 2(a), the ROMs must be of the same dimension and all bases have to be stored. c. Interpolation among the local frequency response functions (FRFs) [13,14]. We will show later in the paper that interpolating Loewner ROMs [15] built from the local FRFs is equivalent to interpolating the local FRFs. Therefore, interpolation among the local FRFs can be seen as a special case of interpolation among the local reduced-order state-space matrices. Furthermore, we propose a further technique to compress the Loewner ROMs to save the storage space. Although this method is intuitive and easy to implement, it suffers from the problem of fixed poles: the positions of the poles do not change with the parameter, but are rather determined by the parameter values used for interpolation.
The goal of the present paper is twofold: • This paper will propose a pole-matching reduced-order modeling method that interpolates linear multiple-input multiple-output (MIMO) ROMs in the frequency domain. Inspired by modal analysis in mechanical engineering [16], the pole-matching method relies completely on analyzing the positions and residues of the poles, rather than trying to recover the state vector of the FOM. We propose to first convert all ROMs to a unified realization, namely the pole-residue realization that stores the positions and residues of poles explicitly. Then, we match the poles according to their positions and residues in order to capture the evolution of the poles in the parametric system. Finally, we interpolate the positions and residues of all matched poles to obtain the parametric ROMs. This method does not require the storage of basis matrices, is capable of interpolating ROMs of different sizes, and works even when a parametric FOM does not exist, e.g., when ROMs are built by a data-driven MOR method or when FOMs at different parameter values result from different discretization methods or different grids. It can also interpolate ROMs of different nature, e.g., interpolating a ROM built by a mathematical MOR method and another ROM obtained from physical reasoning. It is relatively insensitive to the number of parameters and does not assume specific properties of the FOM, e.g., the affinity property required by the reduced basis method [9]. Numerical results show that the pole-matching method is more accurate than the previously proposed methods for our test cases. • The other goal of this paper is to show the importance of the realization of a dynamical system w.r.t. interpolation. For comparison purposes, we will develop another MIMO PMOR method in the frequency domain, namely the interpolation method for Loewner ROMs. We will show that interpolating the ROMs built by the Loewner framework in the original form is equivalent to interpolating the FRFs. Furthermore, we will also discuss the interpolation of Loewner ROMs in the compressed form, which is more efficient w.r.t. computation and storage. Unlike the pole-matching method, which captures the change of positions and residues of the poles, the interpolation of Loewner ROMs builds parametric ROMs with fixed pole positions. Therefore, using different realizations, the parametric ROMs follow different evolu-tion paths. Although both methods have clear physical meanings, the interpolation method for Loewner ROMs follows the path that the real-world system is unlikely to take, while the pole-matching method follows a more natural path and provides accurate parametric ROMs for all our test cases.
Let us emphasize here that the proposed method is more a reduced-order parametric modeling approach than a new PMOR method. Our pole-matching approach assumes the availability of locally valid surrogate models of the FOM. These are assumed to be obtained at feasible sampling points in parameter space, but they can result from various surrogate modeling methods, like • projection-based (or any other computational) MOR methods that compute a nonparametric ROM at the fixed parameter value, • data-driven approaches like the Loewner framework or dynamic mode decomposition [17], • system identification methods, • etc.
We do not even assume that the local surrogates that we interpolate are obtained from the same approach-we can employ a mixture of surrogate models obtained by any of the methods listed above. A particular suitable area of application for our method would be the situation when only an oracle is available, e.g. a running code producing either a state-space model given a fixed parameter or an input-output sequence of time series or frequency-response data. Our approach is fully non-intrusive as it does not require any further knowledge on the system!

Notation
Throughout the paper, ı is the imaginary unit, M T represents the transpose of the matrix M, and M α,T β denotes (M α β ) T .

Background
In this section, we will briefly review two different types of (P)MOR methods: projectionbased methods and the Loewner framework, which is a data-driven (P)MOR method. Though we do not explicitly use any of these methods in our approach, we will frequently refer to them and will also use them for comparison purposes in the numerical examples section. Therefore, we include this brief review for better readability.

Projection-based (P)MOR
This paper focuses on PMOR of state-space systems in the frequency-domain: where E(p), A(p) ∈ R n×n , B(p) ∈ R n×m I , C(p) ∈ R m O ×n , u(s) ∈ R and p ∈ D. A projection-based PMOR method [7,18] first builds two bases Q, U ∈ R n×k (normally, k n), and then approximates X(s, p) ≈ Ux(s, p) (x(s, p) ∈ R k ) in the range of U , and finally forces the residual to be orthogonal to the range of Q to obtain the pROM: where Now we discuss some subcategories under the general framework presented above: • (P)MOR based on projection using local bases These methods build the bases Q and U only using the data computed at a given parameter value, say p 0 . The resulting ROM is valid for p 0 , but if derivative information with respect to p is also included in Q and U , we can obtain a local pROM [8,19,20]. This type of (p)ROMs can be used as the building blocks in our proposed pole-matching PMOR method. When derivative information is included in the bases, we can achieve Hermitian interpolation in our proposed method, which has the potential of reducing the number of needed parameter samples in order to cover a specific region in the parameter space. • PMOR based on projection using global bases. Many PMOR methods build pROMs by projecting the nonparametric ROM onto a global subspace that contains the data obtained at different parameter values of p, namely p 1 , p 2 , …, p j [7,8,11]. Denoting the bases built by a nonparametric MOR method at p i by U i , the global bases U can be obtained by computing the SVD of [U 1 , U 2 , . . . , U j ]. • PMOR based on interpolating local bases These methods interpolate the bases precomputed at different p values, say U i for p i , to compute the bases for the requested parameter value p * [10,12]. A straightforward interpolation normally does not work due to two reasons: 1. The bases U 1 , U 2 , …, U j may be inconsistent. Suppose that the parametric dynamics is well captured by the family of subspaces U(p) := colspan{U (p)}. Further, let U 1 , U 2 , …, U j well represent U(p) at p 1 , p 2 , …, p j . Then, intuitively, the interpolation of these subspaces is meaningful. Nevertheless, the basis matrices U i for these subspaces computed by a MOR method using sampling at p i will generally yield reduced-order models with states living in different coordinate systems as for any orthogonal K ∈ R k×k , we have colspan{U i K } = colspan{U i } = U i . Hence, if we just interpolate the computed bases U i−1 and U i , we consequently interpolate states in different coordinate systems which leads to inconsistencies and poor numerical approximation in state-space. Therefore, in the general case, one should try to compute consistent bases before we conduct interpolation. In [12], it was proposed to compute the bases U i K for p i "as consistent as possible" with the basis U i−1 for p i−1 , by aligning the subspaces via solving the optimization (Procrustes) problem 2. Directly interpolating orthonormal matrices U 1 , …, U j normally does not result in an orthonormal matrix. For example, if U 2 = −U 1 , a direct linear interpolation at p 1 +p 2 2 gives 0, which is apparently not a basis. Therefore, the success of these methods require interpolating on the correct manifold, e.g., on the Grassmann manifold or the Stiefel manifold.
• PMOR based only on (nonparametric) ROMs These methods build the pROM totally based on ROMs, which may be not generated by projection-based (P)MOR methods. Even if they are generated by projection-based (P)MOR methods, it is assumed that the bases U and Q are unknown. For the ROMs built by data-driven methods, e.g., the Loewner framework, we are in this situation.

MOR using the Loewner framework
The objective of MOR using the Loewner framework [15] is to construct a ROM in the frequency domain in the form of using measured samples of the FRF: The Loewner framework assumes that the FRF is only known as a black box: the matrices E, A, B and C are assumed to be unknown but the FRF can be measured at frequency samples. For MIMO systems, each parameter sample is paired with either a left or right tangential direction for the measurement: 1. A right triplet sample: (λ i , r i , w i ). Given a frequency sample λ i and the right tangential direction r i ∈ C m I ×1 , we measure w i = H(λ i )r i . 2. A left triplet sample: (μ j , j , v j ). Given a frequency sample μ j and the left tangential direction j ∈ C 1×m O , we measure v j = j H(μ j ).
Given n R right samples and n L left samples, the Loewner framework computes the Loewner Matrix L ∈ R n L ×n R with and the Shifted Loewner Matrix L σ ∈ R n L ×n R with Using the Loewner matrix and the shifted Loewner matrix, a ROM can be constructed as follows [15].
• The Loewner ROM in the original form can be constructed as: It is highly likely that the matrix pencil (L σ , L) is (numerically) singular, but even in that case, the original Loewner ROM defined in (8) serves as a singular representation of an approximated FRF.
• The Loewner ROM in the compressed form is a concise and regular ROM computed from the Loewner ROM in the original form. First, we compute a rank-revealing SVD in a compressed form: where s is a frequency sample freely chosen from the set {λ i } ∪ {μ j }, and k is the number of dominant singular values chosen for the truncated SVD. Then, a Loewner ROM in the compressed form is constructed as [15]: The Loewner framework has been extended to accommodate parametric systems as a data-driven PMOR method [21]. However, in "Interpolatory PMOR in the Loewner realization" section, we will study another possibility, namely the interpolation of Loewner ROMs.

The pole-matching PMOR method based on the pole-residue realization
In this section, we will first introduce the pole-residue realization for ROMs and develop a pole-matching PMOR method for single-input single-out (SISO) systems, part of which was covered in our conference paper [22]. Then, we will generalize the pole-matching PMOR method to interpolate MIMO ROMs in "The pole-matching method for MIMO systems" section.
The pole-matching method relies exclusively on ROMs at samples where is built for the parameter value p i (i = 1, 2, . . . , n p ). It is not required that the same MOR method is used to build ROMs for all p i values. However, when we do apply a projection-based MOR method to the FOM (1) to obtain ROMs in the form of (2) at each parameter sample p i , it is easy to obtain the ROMs that we need here in the form of (11) as long as E(p i ) is nonsingular by assigning

Remark 1
The assumption of nonsingular E(p i ) is satisfied in many important cases. For example, it always holds when we apply a projection-based MOR method to reduce a system of parametric ordinary differential equations (ODEs) (1) at p i because E(p i ) = Q T E(p i )U , Q and U are both of rank k, and E(p i ) is nonsingular in a system of ODEs. In many situations, e.g., when the model stems from a parametric finite-element model, E will be constant and nonsingular ((projected) mass matrix), and this is the typical class of models we are considering here. For such models, a parameter-dependent E(p) may occur if the geometry of the domain on which the finite-element model is constructed is parameterized. Still, then E(p) will be a mass matrix and will in general be nonsingular, also after (Petrov-)Galerkin projection. Even for differential-algebraic equations with differentiability index one, classical model reduction techniques like balanced truncation [18] or IRKA [23] usually yield reduced-order models with nonsingular E by allowing a nonzero feedthrough-term Du in the output equation; see [24] for details.

The pole-residue realization for SISO systems
Before presenting the MIMO case in "The pole-matching method for MIMO systems" section, we focus on the pole-matching method for SISO systems, i.e., systems in the form of (11) with r I = r O = 1, and for simplicity of notation, we denote B j = B j,1 and C j = C 1,j . In this section, we focus on a single ROM built at the parameter sample p i and omit the index · (i) in the system (11) for simpler notation: For now, we assume that the matrix A is nonsingular and all its eigenvalues are simple: the more complicated cases will be discussed in "Practical considerations" section. For a real eigenvalue λ j and its corresponding eigenvector v j , we have while for the conjugate complex eigenpairs (a j ± ıb j , r j ± ıq j ), the definition A(r j ± ıq j ) = (a j ± ıb j )(r j ± ıq j ) leads to: where for the single real eigenpair (λ j , v j ), while for the conjugate complex eigenpairs (a j ± ıb j , r j ± ıq j ), To facilitate the later more generic discussion for semisimple and defective eigenvalues, we assume j ∈ R m j ×m j and P j ∈ R k×m j in general with m j a possibly larger integer. Then, the complex eigenvalue decomposition is described by the following real matrix Eq. [25]: and it follows the associated similarity transformation Therefore Define where For the real eigenpair (λ j , v j ), C I j and B I j are scalars, with which we define and derive while for the conjugate complex eigenpairs (a j ± ıb j , r j ± ıq j ), we first define , and then derive where we define Therefore, where Definition 1 For the linear system (A, B, C) in (11), its pole-residue realization is defined as ( , B II , C II ) in (25), namely The pole-matching method The pole-residue realization is a natural choice for the interpolation of ROMs because of the following theorem. Proof This is an apparent result from construction since in the pole-residue realization, the positions and residues of the poles are directly stored in the matrix and C II , respectively, for both real eigenvalues and conjugate complex eigenvalues, as is shown in Eqs. (14), (15), (21), (22), (23) and (24).
Remark 2 This theorem shows that when we interpolate the pole-residue realizations ( (i) , B II,(i) , C II,(i) ), we interpolate the positions and residues of the poles, which are values intrinsic to the FRFs themselves. By contrast, interpolating systems of the realization ( (i) , B I,(i) , C I,(i) ), where we only assume that (i) takes the form of (13) and impose no requirement on B I,(i) and C I,(i) , may introduce additional degrees of freedom introduced by realization. For example, assume that ( (i) , B I,(i) , C I,(i) ) (i = 1, 2) are two ROMs of dimension 3 describing the same system at the same parameter value with different realizations: (1)  . These two ROMs have the same FRF due to (22) and the residues for all poles are 16. A critical property of a sound interpolation-based framework is that, when we interpolate two equivalent systems at the same parameter value, the interpolated system should be equivalent to these two systems no matter what interpolation method we use. However, if we interpolate them linearly with the equal weight 0.5, we get (3) = diag{λ 1 , λ 2 , λ 3 }, B I,(3) = [10, 3, 2.5] T , C I,(3) = [2.5, 6, 10], and the residues of the three poles become 25, 18 and 25, respectively, which are all wrong. Therefore, the motivation to define the pole-residue realization is to remove the additional degrees of freedom, which may give rise to spurious results, by "modifying" C I and B I to C II and B II .

Remark 3
The computational cost of the eigendecomposition used in computing the pole-amplitude realization is low since we apply it to a ROM, which is typically of order O(10), in most cases less than 100. Therefore, considering the computing cost alone, any eigendecomposition method can be used. For numerical issues, we refer to "On defective eigenvalues" section because we will discuss general cases there without making the assumption that all eigenvalues of A are simple.
Remark 4 Although we have employed the eigendecomposition to compute the poleamplitude realization like in modal analysis, the eigenmodes of the system may not be well approximated. Nevertheless, a ROM will not accurately capture the dynamics of the FOM unless it captures the dominant eigenmodes well enough, so usually, the dominant eigenmodes are captured quite well in a good ROM. Since our starting point is the ROM rather than the FOM, our full focus is on the input-output behavior of the system. Therefore, the computation of the modes of the system, which depends on all state variables of the FOM, is beyond our consideration. If eigenmodes are to be preserved, then it would be suggested to compute the ROMs at p i using modal truncation [26] so that at p i , the local ROMs contain the exact dominant eigenmodes. By our interpolation approach, the eigenmodes at non-sampling parameters are then approximated by the interpolated poles so that we expect good approximation up to the unavoidable interpolation error.

On the storage and interpolation of ROMs
Besides the state-space representation, the pole-residue realization can be stored and interpolated in a more efficient scheme.
• For a real eigenpair (λ j , v j ), we only need to store two real numbers: λ j and C II j = C I j B I j because B II j ≡ 1 does not need to be stored. • For a complex eigenpair (a j ± ıb j , r j ± ıq j ), we only need to store four real numbers: Therefore, for an order-k ROM in the pole-residue realization, the storage is only 2k: the vector (λ j , C II j ) for a real eigenvalue and (a j , b j , C II j,1 , C II j,2 ) for two conjugate complex eigenvalues.
Assuming that we have n s real eigenvalues and n d pairs of conjugate complex eigenvalues, we store the ROM with two matrices: where each row of D stores a vector of the form (a j , b j , C II j,1 , C II j,2 ) and each row of S stores a vector of the form (λ j , C II j ). Besides the storage efficiency, this storage scheme also has the following advantages: • The order of eigenvalues can be easily rearranged, which provides us great flexibility in the pole-matching process. • An unimportant pole can be easily removed, e.g., using the concept of pole dominance [27]. This can be useful when we want to interpolate ROMs of different orders: for example, when a pole in one ROM cannot be matched to any pole of the other ROM and it is of low dominance, it can be removed in the interpolation process. • It can be easily written back into the state-space representation.
Under the assumption that the matrix A is nonsingular and all its eigenvalues are simple, the pole-matching process to evaluate the ROM at p * / ∈ {p 1 , p 2 , . . . , p n p } works as follows: 1. Given n p ROMs built at p 1 , p 2 , . . . , p n p , we first convert all these ROMs into the pole-residue representation {D (i) , S (i) }. 2. To get the ROM for p * , we first choose an interpolation algorithm and accordingly, the pre-computed ROMs built at p's near p * . For simplicity of presentation, we assume that p 1 and p 2 are chosen for the interpolation. 3. Match the positions and residues of the poles by matching the rows of D (1) and D (2) , and the rows of S (1) and S (2) , respectively. Denote the models after pole matching by D M and S while for a conjugate complex eigenpairs (a (1) where w is a positive real number for weighting. 4. A global "merit function" can also be used, in which case the sum of the local "merit functions" is minimized. Suppose that the two ROMs have the same numbers of real poles and complex poles, respectively. We fix the order of the first ROM and represent the order of the poles in the second ROM after pole-matching by the vector ν, which is a permutation of (1, 2, . . . , n d + n s ). Then we solve the following optimization problem to find ν: Note that all these methods have limitations. The first three methods may result in conflicts in pole-matching, a pole of one ROM is matched to multiple poles of the other ROM.
Since the poles move when the parameters change, it can happen that λ 1 (p 1 ) ≈ λ 2 (p 2 ), λ 1 (p 2 ) ≈ λ 2 (p 1 ), and λ 1 (p 1 ) is far from λ 1 (p 2 ), which we call pole-crossing. If we simply match λ 1 (p 1 ) with λ 2 (p 2 ), and λ 1 (p 2 ) with λ 2 (p 1 ), we lose the true parametric dynamics of the problem. Pole crossing cannot be captured by the first and second methods, and the third and fourth method can also fail. But a trial and error method in pole-matching can always be tried if the engineer is able to tell whether the interpolated ROM is physically sound. When we are capable to build ROM ourselves, e.g., we have access to the FOM, we can build ROMs at more samples to better exploit the parameter space. An offline-online method can be developed to overcome these difficulties, but that is future work. However, the current paper focuses on the cases where the ROMs are given and the following requirements are fulfilled: 1. The given ROMs are accurate enough. Otherwise, we may not be able to match the poles even when the change of parameters is very small. 2. Sufficiently many ROMs are provided to represent the parametric dynamics of the system. To compute the ROM at p * , the poles of the ROMs chosen for interpolation should not change too much. Otherwise, pole matching is difficult due to the lack of data.

Practical considerations
For simplicity of presentation, the discussion above assumed that the matrix A in (11) is real with simple eigenvalues. Now we extend the method to more general cases.

The case of complex A
Assume that A in (12) is a complex matrix with simple eigenvalues. Now in general, the complex eigenvalues occur no longer in conjugate pairs. Therefore, we simply conduct the complex eigenvalue decomposition to diagonalize A. The following computational procedure and storage scheme are the same as in the case of real A with all eigenvalues real. Note that when A is complex, the "residues" stored in C II j are also complex numbers in general. Therefore, we need to interpolate both the positions and the "complex residues" of the poles.

On semisimple eigenvalues
Assume that the dynamical system (11) has a semisimple eigenvalue λ j with multiplicity m j . Then j in (13) is an m j × m j diagonal matrix with all diagonal elements λ j , and C I and B I are row vector and column vector of length m j , respectively. Its corresponding contribution in the sum (20) is As this derivation shows, a simple solution is to define and treat it as if it were a simple pole, i.e., we need only to store (λ j , C II j ) since B II j ≡ 1. However, in reality, the simple strategy does not always work. Consider a parametric system with two poles, one pole with multiplicity 2 for any parameter value and the other pole normally simple. If the two parametric poles coincide at the parameter value p * , the simple procedure above just gives one pole with multiplicity 3, which introduces difficulty in interpolation with ROMs built at other p values, which have two poles. Therefore, we must "separate" the two poles even they happen to be at the same position. This problem will be discussed in more detail in our future work.

On defective eigenvalues
When A has defective eigenvalue(s), the dynamical system (11) cannot be written into the proposed pole-residue realization because the eigenvalue decomposition (16) no longer exists. In this case, A is similar to a Jordan matrix J = P −1 AP, the numerical computation of which is highly unstable [25]. Although in practical computations, defective eigenvalues rarely occur, especially for a ROM, due to numerical noise, a nearly defective eigenvalue can also lead to P having a very large condition number, which may cause numerical instability, e.g., in the computation of B I = P −1 B in (19). Therefore, in practical computations, we always check the condition number of P, which can be computed by, e.g., the MATLAB function "condeig". When it is very large, the algorithm breaks and fails. The solution of this problem is future work.

The pole-matching method for MIMO systems
Now we generalize the pole-matching PMOR method to MIMO systems. When r I > 1 and/or r O > 1 and all eigenvalues of A ∈ R k×k are simple, our derivation in "The poleresidue realization for SISO systems" section holds until (19) with C I ∈ R r O ×k , B I ∈ R k×r I , C I j ∈ R r O ×m j , and B I j ∈ R m j ×r I .

Eigenpair (λ j , v j ) with m j = 1
To study an individual term C I j (sI − j ) −1 B I j in (20), let f denote the index of the first non-zero entry of B I j (which must exist, since otherwise the pole can be removed), and define The contribution of (λ j , v j ) in the weighted sum (20) is Remark 5 As we have discussed in Remark 2, C I j and B I j are not uniquely defined for a given FRF because of the realization freedom. However, the positions and residues of the poles depend only on the FRF. In the case above, the position of the pole is λ i and the residues of the MIMO system are the entries of the matrix Therefore, all entries of C II j are determined by the residues of the poles. Actually, all entries of B II j are also determined by the residues of the poles because where C I j,g is any nonzero entry for C I j . However, B II j,i cannot be used for interpolation purposes. For example, if we interpolate a ROM built for p 1 and another ROM built for p 2 with their weights ω(p) and (1 − ω(p)), respectively, we should compute B II j,i by , which first estimates the residues of the poles by interpolation and then compute B II j,i with these residues, rather than by , which loses the connection with the residues of poles at p. Therefore, it is insufficient to only store B II j for interpolation purposes. We need to store and compute B II j as

The pole-residue realization for MIMO ROMs
For the complex eigenpairs (a j ± ıb j , r j ± ıq j ) with m j = 2, it is difficult to derive where C II j ∈ R r O ×2 and B II j ∈ R 2×r I with all their entries either constants or the residues of the poles. This requires solving a system of nonlinear equations under constraints, which is difficult. Actually, we do not even know whether the solution exists in general. The research of this topic is future work. Therefore, we propose two methods to obtain the pole-residue realization for MIMO ROMs.
The pole-residue realization for MIMO ROMs in the complex form When it is not required to preserve real realization for real systems, a pole-residue realization can be easily computed if all eigenvalues are simple. We substitute the similarity transformation (17) with a true eigendecomposition with diagonal rather than block diagonal. Therefore, we can apply the method developed in "Eigenpair (λ j , v j ) with m j = 1" section to all the eigenvalues to compute the pole-residue realization.
An advantage of this method is that the MIMO structure is strictly preserved: the sizes of , B II and C II equal those of A, B and C, respectively. However, for a real system with conjugate complex eigenpairs, complex numbers are introduced in the pole-residue realization.
The pole-residue realization for MIMO ROMs in the real form To derive the pole-residue realization for MIMO ROMs in the real form, we first consider the SIMO (single-input multiple-output) case.
The SIMO case For a SIMO system (12) with  (17) is the same because for all these SISO systems, the positions of all poles are the same, respectively. Therefore, [i] = . • According to (21) and (24) Due to the property (39), the pole-residue realization for the SIMO system is ( , B II , C II ) with The MIMO case For a MIMO system (12) with C ∈ R r O ×k , B ∈ R k×r I , we denote where B[i] denotes the i-th column of B.
We first compute the pole-residue realization for each SIMO system (A, B[i], C), which we denote by ( , B II , C II {i}). Due to a similar argument as that for (39), we deduce that and B II does not change with i. However, C II {i} does change with i because under different input vectors and output vectors, the residues of each pole are different in general.
Using the method proposed in [28], which reformulates the MIMO systems as parallel connection of split systems, the pole-residue realization of the MIMO system (12) in the real form is = diag , , . . . , , Although this realization preserves the MIMO structure with real algorithms, the dimension of the ROM is multiplied by the number of inputs. In practical computations of FRFs, however, we do not need to formulate (42) explicitly. We compute the FRF column-wise, i.e., to compute the i-th column of the FRF, we only need to formulate ( , B II , C II {i}) for (A, B[i], C).

On the storage and interpolation of MIMO ROMs
In the MIMO case, the FRF H(s) is a matrix function with r O × r I entries. We denote the function for the (j, k)-th entry by H j,k (s). All these functions have the same poles, which need to be stored only once. Normally, they have different residues for the poles, which needs to be stored individually. To conduct ROM interpolation, we conduct pole matching and pole interpolation using similar methods as discussed for the SISO case. A major difference is that now we have more residue information, which we can use to construct the merit function. With the positions and residues of all poles for all FRF entries, an interpolated ROM can be constructed. The previous two sections actually give two procedures that construct a ROM from the positions and residues of all poles, one for a complex realization and the other for a real realization.
For the complex realization discussed in "The pole-residue realization for MIMO ROMs in the complex form" section, we only need to store the diagonal elements of and the matrices C II in (31), along with B II,U j and b II,L j (j = 1, 2, . . . , k) in (35), a total of a total of k × (r O + r I + 2) complex numbers. When we conduct interpolation using this realization after pole-matching, we also need to check that the indices g and f in (35) are the same for each j, respectively.
For the real realization form discussed in "The pole-residue realization for MIMO ROMs in the real form" section, we need to store the positions of all poles and C II defined in (42), a total of k × (r O × r I + 1) real numbers.

Interpolatory PMOR in the Loewner realization
In this section, we propose a method that interpolates ROMs built by the Loewner framework to approximate the dynamical system described by the FRF H(s, p). We will deal with the interpolation of Loewner ROMs in the original form in "Interpolating Loewner ROMs in the original form" section, and the interpolation of Loewner ROMs in the compressed form in "Interpolating Loewner ROMs in the compressed form" section. Both methods rely on Assumption 1.

Assumption 1
We assume that for all sampled values for the parameter p l (l = 1, 2, . . . , n p ), the same frequency samples and right/left tangential directions are used. Therefore, given the frequency shifts μ i and the corresponding left tangential direction i (i = 1, 2, . . . N ω ) a left sample of H(s, p l ) is defined by where v i (p l ) = i H(μ i , p l ).
Similarly, given the frequency shifts λ j and the corresponding right tangential direction r j (j = 1, 2, . . . N ω ), a right sample of H(s, p l ) is defined by Under Assumption 1, the Loewner matrix and the shifted Loewner matrix at p l , which we denote by L(p l ) and L σ (p l ), respectively, are defined by

Interpolating Loewner ROMs in the original form
The following theorem shows the physical meaning of interpolating Loewner ROMs in the original form. and (λ j , r j , w j (p l )), which satisfy Assumption 1. Then, the following two pROMs built for the parameter value p are equal: a. The pROM (E a (p), A a (p), B a (p), C a (p)) obtained by applying an arbitrary interpolation operator of the form to each of E, A, B and C. b. The pROM (E b (p), A b (p), B b (p), C b (p)) built by the Loewner framework in the original form using the "interpolated left/right data" at p, which is obtained by using the interpolation operator (46) on the left/right samples of the FRF: A proof for Theorem 2 is given in Appendix A. Theorem 2 shows that interpolating Loewner ROMs in the original form results in a pROM that is "optimal" from the perspective of the left/right triplet samples. However, this method is practically unsatisfactory as it needs more memory storage than the original left/right triplet samples. Therefore, our ultimate goal is to interpolate Loewner ROMs in the compressed form.

Interpolating Loewner ROMs in the compressed form
To study the interpolation of Loewner ROMs in the compressed form, we first parameterize Eq. (9) by denoting the (truncated-)SVD at the parameter value p l by For a proof of Proposition 1, we refer to Theorem 5.2 in [15]. Therefore, we can compress the ROM in the original representation by ignoring the hardly controllable and observable vectors from X l and Y l , i.e., taking the first k dominant columns of X l and Y l , respectively, and then projecting the state vector to the range of X l,k and the dual state vector to the range of Y l,k . This procedure leads to a Loewner ROM in the compressed form as in (10). Now we propose Algorithm 1 for generating Loewner ROMs in the compressed form that can directly be used for interpolation.
where Y, H ∈ R n×n , X H ∈ R n×n p n , and Y K , H,K and X H,K are obtained by truncated SVD using the first K singular values. 2: Build the global basis X by computing the economy-size SVD of where X, V ∈ R n×n , Y ∈ R n p n×n and Y V,K , V,K and X K are obtained by truncated SVD using the first K singular values. 3: Build the "Compressed" Representations using the global bases: A proof for Theorem 3 is given in Appendix A. According to Theorem 3, if we truncate X an Y to eliminate the hardly controllable and observable subspaces, respectively, all ROMs for l = 1, 2, . . . , n p are accurately approximated by (51).

Remark 6
The Loewner pROM in the compressed form (51) is a good approximation of the Loewner pROM in the original form (46), as long as K is large enough so that X K and Y K capture all dominant components of X and Y . This is because • The pROM (51) is actually obtained by applying the projection method with global bases to the pROM (46). • At each interpolation point p l , the controllability and the observability is captured well by (51) according to Theorem 3, i.e., the Loewner pROM (51) in the compressed form interpolates the Loewner ROMs in the original form well.

Results
In this section, we apply the developed methods to three applications: a microthruster model [19], a "FOM" model [21], and a footbridge model [29]. We compare three methods: the pole-matching method, interpolatory PMOR in the Loewner realization, and the interpolation of ROMs on the nonsingular matrix manifold [12].

PMOR on the microthruster model
In this section, we study the performance of the proposed methods on data-driven ROMs in the frequency domain. The FRFs used to compute the data-driven ROMs are generated by the microthruster model [19]. A schematic diagram of the microthruster is shown in Fig. 1. The microthruster model is of the form (4) with order n = 4257 and has a single parameter: the film coefficient. First in Fig. 2, we show the convergence of nonparametric Loewner ROM in the compressed realization, which is generated with 100 samples for λ and 100 samples for μ. With the increase of the dimension, the FRF of the Loewner ROM becomes closer to that of the FOM. Therefore, the Loewner ROMs are suitable for our study of the interpolation of ROMs.
Then, we apply the three PMOR methods to the microthruster model.
• The manifold method This method interpolates ROMs on the nonsingular matrix manifold [12]. In the numerical tests, we first build ROMs at the parameter samples  Fig. 3. The approximation quality improves as the order of the ROM k increases up to 10. However, when k > 10, the approximation quality becomes unacceptable. Because of its unsatisfactory performance, we will not test this method any further in the next two numerical examples. • Interpolation of ROMs in the Loewner realization (Algorithm 1) First, we interpolate the Loewner ROMs in the original form. Figure 4 shows that this method is much more accurate than the manifold method.
Furthermore, the pROM is much more stable: we have never observed divergence as we increase the order of the pROM. Then, we interpolate the Loewner ROMs in the compressed form built by Algorithm 1. In Fig. 5, we show the numerical results when different shifts s ∈ {λ i,l } ∪ {μ j,l } (defined in (48)) are used. Within each subfigure, the nonparametric ROMs used for interpolation are built for all p l 's using the same shift s specified in the subtitle. It is shown that no matter what shift we use, the resulting ROM is accurate. As a reference, we show in Fig. 6 the numerical results for interpolating the Loewner ROMs in the compressed form using local bases (10) rather than using the global bases in Algorithm 1. The numerical results show that using individual bases rather than global bases in compressing the ROMs, we cannot obtain a pROM with high fidelity by interpolation. In Fig. 7a, we plot the response In Fig. 7b, we show that using a more advanced spline interpolation, higher accuracy can be achieved. • Interpolation of ROMs in the pole-residue realization To apply the interpolation method based on the pole-residue realization, we first study the pole-matching criterion. We use the criterion that the overall differences in pole positions and differences in pole residues are minimal, which agrees with the intuition on the optimal polematching solution. Figure 8 shows the poles of the ROMs for p 22  In this example, the ROMs are of the form (4) with complex E, A, B and C. Since E is nonsingular, we left multiply the system by E −1 , which is of the reduced dimension, to obtain the system of the form (11), for which the pole-residue realization is defined. In this example, the matrix A is complex. To conduct PMOR, we first convert the ROMs (in the form of (11)) at p 1 , p 8 , p 15 , p 22 and p 29 to the pole-residue realization. Then, we conduct linear interpolation among the resulting ROMs. The result is shown in Fig. 7(c). Its accuracy is slightly better than the linear ROM interpolation of the Loewner representation.

Results on the parametric "FOM" model
Now we apply our method to the parametric "FOM" model presented in [21], which is adapted from the nonparametric "FOM" model in [30]: where C = [10, 10, 10, 10, 10, 10, 1, . . . , 1], B = C T , and A(p) = diag(A 1 (p), (1, 2, . . . , 1000), We first use interpolate ROMs in the Loewner realization (Algorithm 1) to obtain a pROM. This example, however, shows the limitation of Algorithm 1 in dealing with peak(s) that moves significantly with the change of parameter(s). As was discussed in [7], Fig. 7 Response surface and the absolute error. a Linear interpolation of the Loewner representation. The overall relative error is 2.4649 × 10 −2 . b Spline interpolation of the Loewner representation. The overall relative error is 6.0671 × 10 −3 . c Linear interpolation of the pole-residue representation. The overall relative error is 2.3485 × 10 −2 when we interpolate FRFs, the positions of poles do not change. This is because the interpolated FRF has a pole at any s, at which any of the individual FRFs has a pole. Therefore, the poles of the interpolated FRF are the union of the poles of the interpolating FRFs. This phenomenon is clearly shown in Fig. 9: with the change of the parameter, the pROM generated by Algorithm 1 does not capture the "moving peak", which is the true dynamics, but rather evolves by waxing and waning of two fixed peaks, which also interpolates the two FRFs but seldom occurs in real applications.

Remark 7
We note the research efforts in avoiding the problem of "fixed peaks" when we interpolate the FRFs. For example, a pair of scaling parameters were introduced in [31]. The method works well when all peaks move in the same direction at a similar rate, which was proved by their numerical tests. In general, however, the method is insufficient to describe the movements of all poles because it actually only introduces one additional degree of freedom.
The ROM interpolation based on the pole-residue realization, on the contrary, is capable of capturing the moving peak because it interpolates the positions of the poles. In this example, we use two MOR methods to build the nonparametric ROMs at the parameter samples: a system identification method as a MOR method (the ssest function in MAT-LAB) [32], and a balanced truncation (BT) method (the balred function in MATLAB). In both cases, the two ROMs for interpolation are built at p = 10 and p = 32.5 with the Fig. 10 The poles of the ROM in the pole-residue realization generated by different MOR methods: the balanced trancation (BT) method and the data-driven method "ssest" in MATLAB dimension k = 10. The positions of the poles of the ROMs at these two parameters are shown in Fig. 10.
In Fig. 11a, b, we interpolate ROMs built by balred and ssest, respectively. The figures show that in this case, balanced truncation achieves better overall accuracy than "ssest". However, a more important observation is that, the errors at the interpolated parameter values are comparable to the errors at the interpolating parameter values. Therefore, the bigger error of the interpolated "ssest" pROM in the pole-residue realization results from the bigger error of the nonparametric "ssest" ROMs used for interpolation, rather than from interpolation itself. So in both cases, ROM interpolation based on the pole-residue realization gives satisfactory results. Now we try to interpolate ROMs of different types. In Fig. 12a, we interpolate a BT ROM built at p = 10 and an ssest ROM built at p = 32.5. From Fig. 10, we can see that the nondominant poles of the FOM are presented by significantly different poles of the ROM in the two different types of methods. However, their interpolation also gives accurate results as Fig. 12a shows. 1 Note that when we interpolate different types of ROMs, we should be particularly careful about pole-matching. If we skip the pole-matching procedure, we will get the result shown in Fig. 12b, which clearly presents the wrong evolution of peaks due to the interpolation of wrong poles.

Results for the footbridge model
In this section, we consider a large-scale footbridge model. The footbridge is located over the Dijle river in Mechelen (Belgium). It is about 31.354 m in length and a tuned mass damper is located in the center. The discretized footbridge model is K 0 + iωC 0 + (k 1 + iωc 1 )K i − ω 2 M 0 X(ω, k 1 , c 1 ) = F, Y (ω, k 1 , c 1 ) = LX(ω, k 1 , c 1 ), where K 0 and M 0 are obtained from a finite element model with 25,962 degrees of freedom, C 0 represents Rayleigh damping, K 1 is a matrix with four non-zero entries that represents the interaction between the tuned mass damper and the footbridge, the input vector F represents a unit excitation at the center span, and the output vector L picks out the vibration at the center span. The model has two parameters, the stiffness of the damper k 1 and the damping coefficient of the damper c 1 . To reduce the model, we apply the Krylov subspace method [33,34] on the first-order equivalent system to obtain ROMs of the form (4) and then, we left multiply the system by E −1 to obtain the system of the form (11). The order of ROMs is set to 10.
In this example, we use four points (k 1 , c 1 ) in the parameter space for interpolation: (10,000 N /m, 20 Ns /m), (20,  The FRFs corresponding to these four points are shown in Fig. 13. Using these four points, we conduct a 2-dimensional linear interpolation (function "interp2" in MATLAB) based on the pole-residue representation to get the ROM for (k 1 , c 1 ) = (15,000 N /m, 35 Ns /m). Figure 14 shows that the interpolated ROM is accurate.

Discussion
Here are some further discussions: • The advantages of the ROM interpolation method based on the pole-residue realization.
• We do not need to know the explicit parametric expression of the system matrices because the parameters only vary in the interpolation of ROMs. • It does not assume the existence of the FOM and works well also with ROMs built by data-driven MOR methods.
• It can even interpolate ROMs built by different MOR methods.
• Its computational cost is relatively insensitive to the number of parameters.
• It can deal with complex parameter dependency, e.g., nonlinear or nonaffine dependence. Since we employ an external interpolation method to handle the  parameter dependency, the proposed method is effective as long as the parameter dependency can be locally captured well by the interpolation method.
• About stability. If we use linear interpolation for the pole-matching PMOR method, the stability is preserved since interpolating the poles in the left-half plane results in poles in the left-half plane. Using other interpolation methods, the stability can in general not be guaranteed. However, we can easily keep track of the interpolated poles in the pole-matching method. A straightforward solution is to use linear interpolation instead when an interpolated pole lies outside the left-half plane. The Loewner PMOR interpolation method always preserves the stability because the FRF at any parameter is a weighted sum of FRFs corresponding to stable systems.

Conclusions
A pole-matching PMOR method that interpolates MIMO linear ROMs in the pole-residue realization was proposed. The method was tested for Loewner-type data-driven ROMs, balanced truncation ROMs, ROMs built by the system identification method ssest, and Krylov-type ROMs. For all these numerical tests, the method gives accurate results. Together with the PMOR method that interpolates MIMO ROMs in the Loewner representation, the importance of realization in ROM interpolation was demonstrated. Similarly, we can prove C a (p) = C b (p). Therefore, the pROM (E a (p), A a (p), B a (p), C a (p)) and the pROM are equal.