## An Adaptive Pole-Matching Method for Interpolating Reduced-Order Models

Yao Yue, Lihong Feng, and Peter Benner,

Abstract -An adaptive parametric reduced-order modeling method based on interpolating poles of reduced-order models is proposed in this paper. To guarantee correct interpolation, a pole-matching process is conducted to determine which poles of two reduced-order models correspond to the same parametric pole. First, the pole-matching in the scenario of parameter perturbation is discussed. It is formulated as a combinatorial optimization problem and solved by a branch and bound algorithm. Then, an adaptive framework is proposed to build repository ROMs at adaptively chosen parameter values, which well represent the parameter domain of interest. To achieve this, we propose techniques including a predictor-corrector strategy and an adaptive refinement strategy, which enable us to use larger steps to explore the parameter domain of interest with good accuracy. The framework also consists of regression as an optional post-processing phase to further reduce the data storage. The advantages over other parametric reduced-order modeling approaches are, e.g., compatibility with any model order reduction method, constant size of the parametric reducedorder model with respect to the number of parameters, and capability to deal with complicated parameter dependency.

FOM is available in the state-space form, which is not always possible in industrial applications. For example, many engineering programs can only compute a FOM from the original partial differential equations at a given parameter value and they may use different meshes for discretization at different parameter values. In addition, they tend to lose efficiency when the number of parameters increases because the dimension of the subspace often increases fast with the number of parameters.

Index Terms -numerical acceleration techniques, reduced order modeling, model-order reduction, frequency domain techniques, spectrum analysis

## I. INTRODUCTION

P ARAMETRIC model order reduction (PMOR) [1] has become a popular tool in the last decades to drastically reduce the computational cost of model-based parametric studies in many fields in engineering, e.g., circuit simulations [2], structural mechanics [3], [4] and chromatography [5].

Among the many methods proposed, projection-based PMOR methods have received the most research attention [6], [7], [1]. Most of these methods first build a global basis by collecting data/derivative information at different parameter values, e.g., snapshots for time-domain systems [6] and (cross-)moments in frequency-domain systems [7]. Then, these methods project the parametric full-order model (FOM) onto the subspace spanned by the basis to compute a parametric reduced-order model (pROM). However, these methods have several disadvantages. First, they assume that a parametric

Y. Yue, L. Feng and P. Benner are with Max Planck Institute for Dynamics of Complex Technical Systems, Magdeburg

Manuscript submitted on July 30, 2019.

This paper is an expanded version from the IEEE MTT-S International Conference on Numerical Electromagnetic and Multiphysics Modeling and Optimization, Boston, MA, USA, May 29-31, 2019.

c 2019 IEEE. Personal use of this material is permitted. Permission from © IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

Recently, there are some research efforts that build pROMs by interpolating local matrices or bases. In [8], it was proposed to interpolate the basis matrices on the geodesic of the Stiefel manifold. However, the costs of both data storage and computation are high and according to our experience in [9], it is difficult to achieve high accuracy with this approach. Another line of research is to conduct PMOR via the interpolation of local matrices. As is shown in [10], directly interpolating the local reduced-order state-space matrices does not work in general since the ROMs are usually in different state-space coordinate systems. Therefore, it was proposed in [10] to compute a global basis from the local bases and to transform all ROMs to a consistent state-space coordinate system defined by the global basis. Another approach was proposed in [11] to apply a congruence transformation to obtain 'consistent bases' in a heuristic sense, which is done by solving a Procrustes optimization problem. It has also been proposed to interpolate the local frequency response functions (FRFs) [12], [13]. However, as we showed in [9], these methods do not work in general. For a comprehensive review of these methods, we refer to [9].

In [14], [9], we proposed a pole-matching based parametric reduced-order modeling method, which builds a pROM by interpolating given nonparametric reduced-order models (ROMs). It interpolates the positions and residues of all poles of the local ROMs, which has a clear physical meaning and removes all additional degrees of freedom introduced by realization. This method only requires these ROMs to be accurate enough, no matter by what methods they are built. We showed that this method can even interpolate a ROM built by a projection-based method and a ROM built by a data-driven method. To match the poles between different ROMs, we proposed to solve a combinatorial optimization problem. However, we did not provide an efficient algorithm to solve this optimization problem: using a brutal-force method quickly becomes too expensive computationally. Furthermore, misleading matching results may occur if the pole distributions of the two ROMs are very different.

In this paper, we improve the efficiency and robustness of

the pole-matching based parametric reduced-order modeling method using the following techniques:

- 1) The branch and bound method for pole-matching. To guarantee correct pole-matching, we propose to use the optimization approach only in the case of pole perturbation, i.e., the corresponding poles of the two ROMs only differ slightly. The proposed branch and bound method is very efficient for this scenario because most enumeration cases are normally branched out thanks to the good initial guess.
- 2) An adaptive framework for parametric reduced-order modeling. This framework serves as an efficient and accurate way to select the parameter values adaptively, at which we build repository ROMs in order to represent the parameter range of interest well. A predictor-corrector technique and an adaptive refinement strategy are proposed to explore the parameter domain of interest using larger-steps with good accuracy. A regression technique is also proposed as a post-processing phase to further reduce the data storage.

This paper is extended from our conference paper [15]. Compared to [15], we propose a new adaptive parametric reduced-order modeling framework, which consists of new techniques like adaptive refinement and regression. The new framework is more accurate and more robust. Our numerical results show that the new method not only improves the relative error from a magnitude of 10 -2 to a magnitude of 10 -4 , but also saves the storage of the pROM. Last but not least, we discuss the branch and bound method in more detail.

This paper is constructed as follows. In Section II, we briefly review the results in [14]. Then in Section III, we discuss matching of perturbed poles and propose the branch and bound method. Based on this, a framework to adaptively build a pROM is proposed in IV. We end the paper with the numerical results in Section V and conclusions in Section VII.

## II. BACKGROUND: THE POLE-MATCHING PMOR METHOD

Assume that we have a series of ROMs and denote the ROM built for the parameter value p i ( i = 1 2 , , . . . , n p ) by ( A ( ) i , B ( ) i , C ( ) i ) :

$$( s I - A ^ { ( i ) } ) x ^ { ( i ) } ( s ) & = B ^ { ( i ) } u ( s ), \\ y ( s ) & = C ^ { ( i ) } x ^ { ( i ) } ( s ), \quad \ \ ( 1 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where A ( ) i ∈ R k × k , B ( ) i ∈ R k × 1 , and C ( ) i ∈ R 1 × k . Assume further that for each i , all eigenvalues of A ( ) i are simple.

## A. The Pole-Residue Realization

Since a ROM has infinitely many equivalent realizations, the simple idea of interpolating A ( ) i , B ( ) i and C ( ) i to obtain a pROM normally does not work [10]. In [14], we proposed to convert all ROMs to a unified pole-residue realization

$$( s I - \Lambda ^ { ( i ) } ) x ^ { ( i ) } ( s ) & = \overline { B } ^ { ( i ) } u ( s ), \quad \ \frac { \varsigma _ { j, 1 } ( s ) } { ( s - } \\ y ( s ) & = \overline { C } ^ { ( i ) } x ^ { ( i ) } ( s ), \quad \ \ ( 2 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where using the notation M α,T β for ( M α T β ) , we define

$$\Lambda ^ { ( i ) } = \text{diag} \left [ \Lambda _ { 1 } ^ { ( i ) }, \Lambda _ { 2 } ^ { ( i ) }, \dots, \Lambda _ { m } ^ { ( i ) } \right ], \text{ \quad \ \ } ( 3 )$$

$$\overline { B } ^ { ( i ) } = \left [ \overline { B } _ { 1 } ^ { ( i ), T }, \overline { B } _ { 2 } ^ { ( i ), T }, \dots, \overline { B } _ { m } ^ { ( i ), T } \right ] ^ { T }, \quad \ \ ( 4 )$$

$$\overline { C } ^ { ( i ) } = \left [ \overline { C } _ { 1 } ^ { ( i ) }, \overline { C } _ { 2 } ^ { ( i ) }, \dots, \overline { C } _ { m } ^ { ( i ) } \right ].$$

The pole-residue realization is closely related to the rational form of the transfer function H ( ) i ( s ) = y s ( ) u s ( ) of (1):

$$\int _ { \substack { n + \text{or} \\ \text{tively}, \\ \text{present} } } H ^ { ( i ) } ( s ) = \sum _ { j = 1 } ^ { n _ { s } } \frac { \overline { C } _ { j } ^ { ( i ) } } { s - \lambda _ { j } ^ { ( i ) } } + \sum _ { j = 1 } ^ { n _ { d } } \frac { \overline { C } _ { j, 1 } ^ { ( i ) } ( s - a _ { j } ^ { ( i ) } ) - \overline { C } _ { j, 2 } ^ { ( i ) } b _ { j } ^ { ( i ) } } { ( s - a _ { j } ^ { ( i ) } ) ^ { 2 } + ( b _ { j } ^ { ( i ) } ) ^ { 2 } }, \, ( 6 )$$

where for a real eigenvalue λ ( ) i j , we have Λ ( ) i j = λ ( ) i j , B ( ) i j = 1 and C ( ) i j ∈ R , which corresponds to the term

$$\frac { \overline { C } _ { j } ^ { ( i ) } } { s - \lambda _ { j } ^ { ( i ) } } = \overline { C } _ { j } ^ { ( i ) } ( s - \lambda _ { j } ^ { ( i ) } ) ^ { - 1 } \overline { B } _ { j } ^ { ( i ) },$$

and for a pair of conjugate complex eigenvalues a j ± ıb j , we have Λ ( ) i j ∈ R 2 × 2 , B ( ) i j ∈ R 2 × 1 , C ( ) i j ∈ R 1 × 2 , which corresponds to the term

$$\begin{array} { c c c } \i i & \quad \cdot \\ \mathfrak { e } & \frac { \overline { C } _ { j, 1 } ^ { ( i ) } ( s - a _ { j } ^ { ( i ) } ) - \overline { C } _ { j, 2 } ^ { ( i ) } b _ { j } ^ { ( i ) } } { ( s - a _ { j } ^ { ( i ) } ) ^ { 2 } + ( b _ { j } ^ { ( i ) } ) ^ { 2 } } \\ \jmath & = \left [ \overline { C } _ { j, 1 } ^ { ( i ) }, \, \overline { C } _ { j, 2 } ^ { ( i ) } \right ] \left ( s I - \left [ \begin{array} { c c c } a _ { j } ^ { ( i ) } & b _ { j } ^ { ( i ) } \\ - b _ { j } ^ { ( i ) } & a _ { j } ^ { ( i ) } \end{array} \right ] \right ) ^ { - 1 } \left [ \begin{array} { c c c } 1 \\ 0 \end{array} \right ] \\ \end{array}$$

$$\stackrel { \triangle } { = } \overline { C } _ { j } ^ { ( i ) } ( s - \Lambda _ { j } ^ { ( i ) } ) ^ { - 1 } \overline { B } _ { j } ^ { ( i ) }.$$

Therefore, the positions of the poles are stored in Λ ( ) i and the residues of poles are stored in C ( ) i . The matrix B ( ) i contains no useful information: B ( ) i j ≡ 1 for a real pole and B ( ) i j ≡ [1 , 0] T for a pair of conjugate complex poles.

For more details, we refer to [14], [9].

## B. Efficient Storage

Assuming that A ( ) i has n s real eigenvalues and n d pairs of conjugate complex eigenvalues, we can store the ROM in the pole-residue realization with two matrices [14]:

$$D _ { i } \in \mathbb { C } ^ { n _ { d } \times 4 } \text{ \ and \ } S _ { i } \in \mathbb { C } ^ { n _ { s } \times 2 }. \text{ \quad \ \ } ( 8 )$$

The j -th row of S i is ( λ ( ) i j , C ( ) i j ) , which corresponds to the real pole C ( ) i j s -λ ( ) i j , while the j -th row of D i is ( a ( ) i j , b ( ) i j , C ( ) i j, 1 , C ( ) i j, 2 ) , which corresponds to a pair of conjugate complex poles C ( ) i j, 1 ( s -a ( ) i j ) -C ( ) i j, 2 b ( ) i j ( s -a ( ) i j ) 2 +( b ( ) i j ) 2 . The order of the rows in D i and S i depends on their order in Λ ( ) i : the order of rows in S i ( D i ) is the order of real (complex conjugate) poles in Λ ( ) i .

## III. MATCHING OF PERTURBED POLES: A BRANCH AND BOUND METHOD

To interpolate two sufficiently accurate ROMs, say ( D ,S 1 1 ) and ( D ,S 2 2 ) , to obtain a pROM, we first determine which poles in ( D ,S 1 1 ) and ( D ,S 2 2 ) correspond to the same parametric pole. This is what we call pole-matching. The order of the poles are influenced by: 1) the MOR algorithm; 2) the eigenvalue decomposition algorithm [16] that we use to compute the pole-residue realization (2). In this section, we discuss matching of perturbed poles under Assumption III.1. The matching of distant poles will be the topic of Section IV.

Assumption III.1. The assumptions for the branch and bound method.

- 1) The sizes of D 1 and S 1 equal those of D 2 and S 2 , respectively;
- 2) All poles of the ROM ( D ,S , 2 2 ) are perturbations of poles of ROM ( D ,S 1 1 ) , but for a pair of perturbed poles, their order in the D S / vectors may be different. More precisely, r ( ROM 1 , ROM 2 ) is small, with r defined later in (26) .

Under Assumption III.1, pole-matching can be done by solving an optimization problem, i.e., find an optimal way to rearrange the rows of ( D ,S 2 2 ) so that it is the closest to ( D ,S 1 1 ) in some sense.

## A. Matching of Perturbed Poles via Optimization

We first define two mapping vectors: v d = [ v d, 1 , v d, 2 , . . . , v d,n d ] for conjugate complex poles and v s = [ v s, 1 , v s, 2 , . . . , v s,n s ] for real poles, which satisfy

$$\{ v _ { d, 1 } \} \cup \{ v _ { d, 2 } \} \cup \cdots \cup \{ v _ { d, n _ { d } } \} = \{ 1, 2, \dots, n _ { d } \}, \quad ( 9 ) \ \text{There}$$

and

$$\{ v _ { s, 1 } \} \cup \{ v _ { s, 2 } \} \cup \dots \cup \{ v _ { s, n _ { s } } \} = \{ 1, 2, \dots, n _ { s } \}. \quad ( 1 0 ) \quad _ { n } ( 1$$

Then we define the mapping M , which rearranges the rows of D (or S ) according to the index vector v d (or v s ):

$$\mathcal { M } ( v _ { d } ; D ) = [ \underset { \begin{smallmatrix} \sim \\ \sim \\ \sim \end{smallmatrix} \pmod { \Phi } ^ { T }, d _ { v _ { d, 2 } } ^ { T }, \dots, d _ { v _ { d, n _ { d } } } ^ { T } ] ^ { T }, \quad ( 1 1 ) \quad \ t h o m \, \alpha$$

$$\mathcal { M } ( v _ { s } ; S ) = [ s _ { v _ { s, 1 } } ^ { T }, s _ { v _ { s, 2 } } ^ { T }, \dots, s _ { v _ { s, n _ { s } } } ^ { T } ] ^ { T }, \quad ( 1 2 ) ^ { \quad t h e n \, a }$$

where D = [ d T 1 , d T 2 , . . . , d T n d ] T and S = [ s T 1 , s T 2 , . . . , s T n s ] T .

In the case of pole perturbation, we achieve pole-matching by solving the following optimization problem

$$\min _ { v _ { d }, v _ { s } } f ( v _ { d }, v _ { s } ) = & \left \| D _ { 1 } W _ { d } - \mathcal { M } ( v _ { d }, D _ { 2 } ) W _ { d } \right \| _ { F } ^ { 2 } + \quad \text{ using } ( 2 1 \\ & \left \| S _ { 1 } W _ { s } - \mathcal { M } ( v _ { s }, S _ { 2 } ) W _ { s } \right \| _ { F } ^ { 2 }, \quad ( 1 3 ) \quad v _ { d, i } ^ { - }$$

where W d = diag { w , w , w p p r , w r } and W s = diag { w , w p r } are weighting matrices. The scalars w p and w r are the weights for the positions and residues of the poles, respectively. Now we give the motivation of the optimization problem. Assume that all poles are simple and the perturbation is sufficiently small. If all poles are correctly matched, the objective function is very small; otherwise, the objective function must be much larger.

## B. A Branch and Bound Algorithm

To solve the combinatorial optimization problem (13), the exhaustive search method requires n d ! × n s ! enumerations, which is computationally feasible only for very small ROMs. In this section, we discuss efficient solution of (13).

First, we observe that the optimization problem (13) can be decoupled into two independent optimization problems:

$$\min _ { v _ { d } } f _ { d } ( v _ { d } ) = \| D _ { 1 } W _ { d } - \mathcal { M } ( v _ { d }, D _ { 2 } ) W _ { d } \| _ { F } ^ { 2 } \quad ( 1 4 )$$

and

$$\min _ { v _ { s } } f _ { s } ( v _ { s } ) = \| S _ { 1 } W _ { s } - \mathcal { M } ( v _ { s }, S _ { 2 } ) W _ { s } \| _ { F } ^ { 2 } \,. \quad ( 1 5 )$$

Their solutions v ∗ d and v ∗ s solve the original optimization problem (13):

$$\overset { \prime } { e s } _ { s i v } \quad \min _ { v _ { d }, v _ { s } } f ( v _ { d }, v _ { s } ) = \min _ { v _ { d } } f _ { d } ( v _ { d } ) + \min _ { v _ { s } } f _ { s } ( v _ { s } ), \quad \quad ( 1 6 )$$

$$\stackrel { \delta, \text{ ther} } { \L c i s e l y }, \quad \arg \min _ { v _ { d }, v _ { s } } f ( v _ { d }, v _ { s } ) = \left ( \arg \min _ { v _ { d } } f _ { d } ( v _ { d } ), \, \arg \min _ { v _ { s } } f _ { s } ( v _ { s } ) \right ). \ ( 1 7 )$$

Using this method alone, we reduce the computational cost of (13) from n d ! × n s ! enumerations to n d !+ n s ! enumerations. In the following discussions, we focus on the solution of (14). The optimization problem (15) can be solved similarly.

The second technique that we use to further reduce the computational cost of the combinatorial optimization problem (13) is the branch and bound method [17]. It can drastically reduce this computational cost in the scenarios that are the most useful in pole-matching, i.e., Scenario 1 and Scenario 2 on page 5. The efficiency of the branch and bound method for pole matching relies on the following theorem.

glyph[negationslash]

Theorem III.1. Assume that 1 ≤ i, j ≤ n d are two given indices ( i = j ). If there exist two mapping vectors v (1) d and v (2) d satisfying (9) and

$$v _ { d, i } ^ { ( 1 ) } = v _ { d, j } ^ { ( 2 ) }, \ \ v _ { d, j } ^ { ( 1 ) } = v _ { d, i } ^ { ( 2 ) },$$

glyph[negationslash]

$$v _ { d, k } ^ { ( 1 ) } = v _ { d, k } ^ { ( 2 ) }, \ \ ( \forall 1 \leq k \leq n _ { d }, \, k \neq i \text{ and } k \neq j ) \quad ( 1 9 )$$

glyph[negationslash]

$$f _ { d } \left ( v _ { d } ^ { ( 2 ) } \right ) - f _ { d } \left ( v _ { d } ^ { ( 1 ) } \right ) = c > 0,$$

then any mapping vector v d satisfying (9) and

$$v _ { d, i } = v _ { d, i } ^ { ( 2 ) } \ \ a n d \ \ v _ { d, j } = v _ { d, j } ^ { ( 2 ) } \ \quad \quad ( 2 1 )$$

cannot be the optimal solution of (14) .

Proof. Assume v + d is the optimal solution of (14) satisfying (21). We construct v -d as glyph[negationslash]

$$( 1 3 ) \quad & v _ { d, i } ^ { - } = v _ { d, j } ^ { + } = v _ { d, i } ^ { ( 1 ) }, \quad v _ { d, j } ^ { - } = v _ { d, i } ^ { + } = v _ { d, j } ^ { ( 1 ) }, \\ w _ { r } \} & v _ { d, k } ^ { - } = v _ { d, k } ^ { + }, \quad & ( \forall 1 \leq k \leq n _ { d }, \, k \neq i \text{ and } k \neq j ).$$

glyph[negationslash]

Then,

$$\int _ { W } ^ { \infty } \sum _ { \substack { f _ { d } ( v _ { d } ^ { + } ) - f _ { d } ( v _ { d } ^ { - } ) } } & = \left \| D _ { 1 } W _ { d } - \mathcal { M } ( v _ { d } ^ { + }, D _ { 2 } ) W _ { d } \right \| _ { F } ^ { 2 } - \\ \text{cint} \\ \text{: much} } & = \left \| \left [ \begin{array} { c } D _ { 1, i } \\ D _ { 1, j } \end{array} \right ] W _ { d } - \left [ \begin{array} { c } D _ { 2, v _ { d, i } ^ { + } } \\ D _ { 2, v _ { d, j } ^ { + } } \end{array} \right ] W _ { d } \right \| _ { F } -$$

Therefore,

$$f _ { d } ( v _ { d } ^ { - } ) = f _ { d } ( v _ { d } ^ { + } ) - c < f _ { d } ( v _ { d } ^ { + } ),$$

which contradicts the assumption that v + d is the optimal solution of min v d f d ( v d ) defined in (14).

Based on Theorem III.1, a branch and bound algorithm is proposed with two motivations to save the computational cost:

- 1) Start at a hopefully good initial guess v (0) d = (1 2 , , . . . , n d ) , which is computationally free.
- 2) Try to decrease f d by swapping any two entries of the current v d and use the result from Theorem III.1 to branch out the swaps that we are certain to be unable to further decrease f d .

Remark III.1. Remarks on the two motivations.

On Motivation 1 . Although v (0) d = (1 2 , , . . . , n d ) is simple, it is actually often a quite good initial guess under the following assumption:

Assumption III.2. Conditions in favor of a good initial guess v (0) d = (1 2 , , . . . , n d ) .

- a) We use the same MOR method to compute the ROMs;
- b) The MOR algorithm normally preserves the order of parametric eigenvalues when ‖ p 1 -p 2 ‖ is small;
- c) Both ROMs are accurate enough and ‖ p 1 -p 2 ‖ is small enough.

Under this assumption, it is reasonable to presume that when the parameter is perturbed, 1) The order of the parametric eigenvalues is normally preserved. 2) Even when the order changes, normally only a small number of indices change at the same time.

On Motivation 2 . Since we start from a 'hopefully good' initial guess, it is quite likely that we fail to decrease f d by swapping entries of v d . This case is actually ideal for us because the failed swap can be branched out and we never need to try it again according to Theorem III.1. If the initial guess happens to be the optimal solution, we actually just need to branch out all possibilities, which only needs n d ( n d -1) enumerations instead of n d ! enumerations.

$$& \left \| \left [ \begin{array} { c } D _ { 1, i } \\ D _ { 1, j } \end{array} \right ] W _ { d } - \left [ \begin{array} { c } D _ { 2, v _ { d, i } ^ { - } } \\ D _ { 2, v _ { d, j } ^ { - } } \end{array} \right ] W _ { d } \right \| _ { F } \quad F \ c \\ = & \left \| \left [ \begin{array} { c } D _ { 1, i } \\ D _ { 1, j } \end{array} \right ] W _ { d } - \left [ \begin{array} { c } D _ { 2, v _ { d, i } ^ { ( 2 ) } } \\ D _ { 2, v _ { d, j } ^ { ( 2 ) } } \end{array} \right ] W _ { d } \right \| _ { F } - \quad 1 \\ & \left \| \left [ \begin{array} { c } D _ { 1, i } \\ D _ { 1, j } \end{array} \right ] W _ { d } - \left [ \begin{array} { c } D _ { 2, v _ { d, i } ^ { ( 1 ) } } \\ D _ { 2, v _ { d, j } ^ { ( 2 ) } } \end{array} \right ] W _ { d } \right \| _ { F } \quad \text{Mo.} \\ = & \left \| D _ { 1 } W _ { d } - \mathcal { M } ( v _ { d } ^ { ( 2 ) }, D _ { 2 } ) W _ { d } \right \| _ { F } ^ { 2 } - \\ & \left \| D _ { 1 } W _ { d } - \mathcal { M } ( v _ { d } ^ { ( 1 ) }, D _ { 2 } ) W _ { d } \right \| _ { F } ^ { 2 } \\ = & f _ { d } \left ( v _ { d } ^ { ( 2 ) } \right ) - f _ { d } \left ( v _ { d } ^ { ( 1 ) } \right ) = c > 0. \\ & \text{def} \quad \text{has}$$

To keep track of the branched out swaps, we use the matrix F defined as follows:

$$F _ { F } ^ { - 1 } \begin{cases} - \quad F ( i, j ) = \begin{cases} 0, & \text{when any $v_{d}$ with $v_{d,i}$ = $v_{d,j}^{(0)}$ and} \\ & v _ { d, j } = v _ { d, i } ^ { ( 0 ) } $ has been branched out; \\ 1, & \text{otherwise.} \end{cases} \\ F _ { F } ^ { - 1 } \begin{cases} 0, & \text{when any $v_{d}$ with $v_{d,i}$ = $v_{d,j}^{(0)}$ and} \\ & v _ { d, j } = v _ { d, i } ^ { ( 0 ) } $ has been branched out; \\ 1, & \text{otherwise.} \end{cases}$$

More specifically, we initially set

$$F = J _ { n _ { d } } - I _ { n _ { d } },$$

where J n d , I n d ∈ R n d × n d are the matrix of ones and the identity matrix, respectively. This means that initially, no swap has been branched out: only swapping an index with itself is forbidden. Note that when we swap the i -th and j -th entries of v d , we actually swap the v d,i -th and v d,j -th rows of D 2 , or equivalently, swap the v d,i -th and v d,j -th entries of v (0) d .

Algorithm III.1. A Branch and Bound Method to Solve the Pole-Matching Optimization Problem (14)

glyph[negationslash]

1: Initialize v d = (1 2 , , . . . , n d ) , the current lowest objective value b l = f d ( v d ) , and F = J n d - I n d . 2: while ‖ F ‖ F = 0 do 3: for i = 1 to n d do 4: for j = 1 to n d do 5: if F v ( d,i , v d,j ) = 0 then 6: continue (go to line 4) 7: else 8: Set v ∗ d = v d . 9: Swap the i -th and the j -th entry of v ∗ d . 10: if f ( v ∗ d , v (0) s ) ≥ b l then 11: Update F v ( d,i , v d,j ) = 0 . 12: else 13: Set v d = v ∗ d and b l = f d ( v ∗ d ) , and break both for loops (go to line 2). 14: end if 15: end if 16: end for 17: end for 18: end while

Now we discuss the definition (22). Assume that we are at v C d and we have failed to increase f d by swapping its i -th entry with its j -th entry. Denote the swapped vector by v S d . Then v C d and v S d satisfy the conditions (18), (19) and (20). Therefore, according to Theorem III.1, any v d with v d,i = v S d,i = v C d,j and v d,j = v S d,j = v C d,i cannot be the solution of the optimization problem (14). Therefore, we set F v ( C d,i , v C d,j ) = 0 to forbid this swap, i.e., if the k -th and the l -th entries of any v d satisfy

$$v _ { d, k } = v _ { d, i } ^ { C }, \ \ v _ { d, l } = v _ { d, j } ^ { C }, \ \ \ \ ( 2 4 )$$

we do not need to try the swap of the k -th and l -th entries of v d since it cannot decrease the objective f d . Note that the matrix F is normally non-symmetric because if F v ( d,j , v d,i ) = 0 , we forbid swapping the k -th and the l -th entries if v d,k = v C d,j and v d,l = v C d,i , which is actually a different swap.

As for implementation, we can also use the more compact 'potentially feasible tables 1 ' F 1 , . . . , F v d as we did in [14] instead of using the matrix F , especially when n d is large. The potentially feasible table is defined as

$$\mathcal { F } _ { i } = \{ j | F ( i, j ) = 1 \}. \quad \quad ( 2 5 ) ^ { \ \mathcal { \lambda } _ { 1 } },$$

So when a swap is branched out, it is simply removed from the potentially feasible table.

Now we discuss several scenarios in pole matching.

Scenario 1 . All poles naturally match: the most optimistic scenario. This case occurs quite often under Assumption III.2 as we have discussed there. Using the branch and bound algorithm, we just need to verify that the initial v (0) d is optimal and we only need to enumerate the v d vector n d ( n d -1) times instead of n d ! times.

Scenario 2 . All poles do not naturally match, but only a few swapping operations, say, m times, are required to obtain the optimal solution. In general, this scenario is unavoidable. e.g., when the dominances 2 of two poles cross. The computational cost in this scenario is bounded by mn d ( n d -1) enumerations. The actual computational cost is normally much lower because many swaps have been branched out at an early stage.

Scenario 3 . Almost all poles between ( D ,S 1 1 ) and ( D ,S 2 2 ) do not naturally match. In this scenario, the computational cost of optimization is high. Actually, if both ROMs are generated by the same MOR method with some 'algorithmic continuity', it is not so meaningful to conduct optimization, whose solution can be misleading. In such cases, it is highly possible that the 'perturbation' is too large, except for some rare cases, e.g., the dominances of many pair of poles do cross at the same parameter value. Therefore, we avoid solving the optimization problem in this case, which is not only computationally expensive, but also useless. This is a motivation of the predictor-corrector method that will be presented in Section IV-A.

## IV. THE ADAPTIVE FRAMEWORK TO BUILD THE PARAMETRIC ROM

The goal of this section is to adaptively build a series of repository ROMs, which will be used as interpolates for constructing the pROM. In this section, we confine our discussion to parametric systems with a single parameter p ∈ R . We denote these repository ROMs, whose poles have already been matched to those of all other repository ROMs, by ROM , ROM , ..., ROM 1 2 N built for the parameter values p 1 &lt; p 2 &lt; . . . &lt; p N , respectively. In this paper, the overline symbol · indicates the name or variables of a repository ROM, e.g., ROM with its pole-residue realization i ( D , S i i ) .

1 We use the term 'potentially feasible' because the swaps recorded in this table might not be feasible, but we have not validated their feasibility yet.

2 In this paper, we use the term 'pole dominance' to denote the importance of a pole on the system output. We use it as a general term for all kinds of MOR methods. This term is borrowed from a specific case: the dominant pole algorithm [18], [4], which defines the dominance of a pole of the form

$$\text{$p$ among $[\iota],[^{\iota}],$ when enimes one common $i$} \\ \frac { \overline { C } _ { j } ^ { ( i ) } } { s - \lambda _ { j } ^ { ( i ) } } \text{ by $\frac{|\overline { C } _ { j } ^ { ( i ) } | } { |\Re \{ \lambda _ { j } ^ { ( i ) } \} | }$ for this specific case.}$$

For ease of notation, we introduce a mapping for polematching P D 1 ( D 2 ) ( P S 1 ( S 2 ) ). When we match D 2 to D 1 , we denote the matching result of Algorithm III.1 by P D 1 ( D 2 ) . Using this notation, we define the distance between ROM 1 ( D ,S 1 1 ) and ROM 2 ( D ,S 2 2 ) as

$$r ( \text{ROM} _ { 1 }, \text{ROM} _ { 2 } ) = & \left \| D _ { 1 } W _ { d } - \mathcal { P } _ { D _ { 1 } } ( D _ { 2 } ) W _ { d } \right \| _ { F } + \\ & \left \| S _ { 1 } W _ { s } - \mathcal { P } _ { S _ { 1 } } ( S _ { 2 } ) W _ { s } \right \| _ { F } \quad ( 2 6 )$$

and the relative error by

$$\Phi _ { \Phi } \Phi _ { 1 }, \text{ROM} _ { 2 } ) = & \frac { \| D _ { 1 } W _ { d } - \mathcal { P } _ { D _ { 1 } } ( D _ { 2 } ) W _ { d } \| _ { F } } { \| D _ { 1 } \| _ { F } } + \\ & \frac { \| S _ { 1 } W _ { s } - \mathcal { P } _ { S _ { 1 } } ( S _ { 2 } ) W _ { s } \| _ { F } } { \| S _ { 1 } \| _ { F } }. \quad ( 2 7 ) \\ \Phi _ { \Phi } \Phi _ { 1 } \Phi _ { 2 } \Phi _ { 2 } \Phi _ { 3 } \Phi _ { 4 } \Phi _ { 5 } \Phi _ { 6 } \Phi _ { 7 } \Phi _ { 8 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phicon \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi _ { 9 } \Phi }$$

For convenience, we denote the pole-matched ROM ( P D 1 ( D 2 ) , P S 1 ( S 2 )) by P ROM 1 ( ROM 2 ) .

## A. The Predictor-Corrector Strategy

The two motivations of the predictor-corrector strategy are:

Objective 1 . We want to design a framework, in which Scenario 1 occurs most often, Scenario 2 occurs sometimes, and Scenario 3 occurs rarely.

Objective 2 . However, we do not want to achieve Objective 1 simply by using small steps in p , which is inefficient.

In order to explore the parameter space adaptively, we make the following assumptions on the ROMs.

Assumption IV.1. At any feasible parameter value p , a ROM of high accuracy can be built at request.

Assumption IV.2. The changes in positions and residues of all poles are small when the change of p is sufficiently small.

Suppose that we have already computed the repository ROMs: ROM , ROM , ..., ROM 1 2 i and we want to expand the series by adding new repository ROMs. At a new parameter value p i +1 ( p i +1 &gt; p i ), we first compute a candidate ROM cand i +1 and try to match its poles to the repository ROMs. Instead of applying Algorithm III.1 directly to ROM i and ROM cand , we i +1 propose to use a predictor-corrector strategy.

We first use the positions and residues of poles of ROM 1 , ROM , ..., ROM 2 i to calculate a 'predicted ROM' ROM pred i +1 for p i +1 . For prediction, we can use extrapolation, e.g., linear, polynomial, spline, etc., on the positions and residues of the poles. If the calculated poles exhibit some noise, we can also use regression. Normally, the poles of ROM pred i +1 are closer to the poles of ROM cand i +1 than those of ROM . To be safe, we i check whether r ( ROM i , ROM cand i +1 ) &gt; r ( ROM pred i +1 , ROM cand i +1 ) holds. If it holds, we set ROM i +1 = P ROM pred i +1 ( ROM cand i +1 ) ; otherwise, we set ROM i +1 = P ROM i ( ROM cand i +1 ) .

Note that although we have 'accepted' the current ROM i +1 as a repository ROM, it is not guaranteed yet that the interpolation will give accurate pMOR on the interval ( p , p i i +1 ) . Actually, the pole-matching is not necessarily correct because Assumption III.1 may not be satisfied. However, even if the poles are wrongly matched, we store it for further use because

we do not want to waste the computational effort with which we built ROM i +1 . Therefore, we introduce an index i h to denote the last index that we are sure to have 'high fidelity', i.e., it is correctly matched and close enough to ROM i h -1 to produce good interpolation results. Here we assign i h = i , which means that the fidelity of ROM i +1 still needs to be checked. More details will be discussed in the next section.

## B. An Adaptive Refinement Strategy

In the previous section, we have discussed the efficient generation and pole-matching of ROMs. However, this alone does not guarantee good accuracy of the pROM: it is possible that the poles are matched correctly, but the interpolation result does not capture the parametric evolution of the poles well. For example, for a system that has only one pole, although the pole is naturally 'matched', we still need a fine grid of p to capture the dynamics of the system when the dynamics is complex.

After we have computed ROM i +1 using the predictorcorrector strategy, we need to assess whether good accuracy can be achieved in the interval ( p , p i i +1 ) using interpolation. To achieve this, we choose the test point p t = p i + p i +1 2 and build two ROMs at p t :

1)

ROM

I,t computed by interpolating the repository ROMs.

- 2) ROM T,t computed by a MOR algorithm.

Since under Assumption IV.1, ROM T,t is highly accurate, we regards ROM I,t as accurate only when

$$e ( \text{ROM} ^ { I, t }, \text{ROM} ^ { T, t } ) < \tau _ { e }, \text{ \quad \ \ } ( 2 8 ) \text{ \ \ } \stackrel { \cdot \cdot } { 7 }.$$

where τ e is a small real number, e.g., 0.001. The value of τ e should take the error of the MOR algorithm into consideration: τ e should be sufficiently larger than an estimated error of the MOR algorithm. Otherwise, the relationship (28) may never be satisfied despite of refinement of the grid. In practical computations, we can actually change τ e adaptively according to an error bound of ROM T,t rather than using the same τ e for all iterations.

If (28) is met, we trust the fidelity of ROM i +1 and use the predictor-corrector for the next step. Otherwise, we need more ROMs to represent the interval ( p , p i i +1 ) . A natural choice is ROM T,t since it has already been built. Therefore, we assign

$$\overline { R O M } _ { i + 2 } \leftarrow \overline { R O M } _ { i + 1 }, \ \overline { R O M } _ { i + 1 } \leftarrow \mathcal { P } _ { \overline { R O M } _ { i } } ( R O M ^ { T, t } ). \ ( 2 9 ) \ \stackrel { \cdots } { 1 9 } ; \ \overline { r }$$

Then, we check the fidelity of the new ROM i +1 using the procedure described above:

- · If (28) is violated, we do a further refinement for the interval ( p , p i i +1 ) .
- · Otherwise, we accept that ROM i +1 has high fidelity, update i h = i +1 and

$$\overline { R O M } _ { i + 2 } \leftarrow \mathcal { P } _ { \overline { R O M } _ { i + 1 } } ( \overline { R O M } _ { i + 2 } ) \quad \ ( 3 0 ) \ \frac { 2 6 \colon } { \gamma _ { 7 } }.$$

because p i +1 (the former p t ) is closer to p i +2 (the former p i +1 ) and the pole-matching result is more trustworthy. Next, we check whether ROM i +2 has high fidelity:

- -If it has high fidelity, we set i h = i +2 and end this local refinement.
- -Otherwise, we do a further refinement for the interval ( p i +1 , p i +2 ) using the procedure above.

For the whole procedure, we refer to Algorithm IV.1.

## C. Compact Parametric ROM by Regression

After we have finished the construction of ROM 1 , ROM , 2 . . . , ROM N for the whole interval of interest, we can further save data storage by using regression. In this paper, we do regression for each entry of D j,k i ( ) (the ( j, k ) -th entry of D i built at p i ) and S i ( j, k ) to obtain the parametric forms D p j, k ( ; ) and S p j, k ( ; ) . Using a polynomial of degree q for each entry, we reduce the storage cost from N (4 n d +2 n s ) to ( q +1)(4 n d +2 n s ) since we only need to store the coefficients of the polynomials. This reduction is usually drastic since a low-order polynomial is normally enough when the system dynamics is not too complex.

We summarize all the techniques discussed in Algorithm IV.1 with more details.

## Algorithm IV.1. An Adaptive Pole-Matching PMOR method

- 1: Input : the parameter range of interest [ p L , p U ] .
- 2: Initialization : initialize step length u 0 , error tolerance τ e , order of regression polynomials q , i = 2 , p 1 = p L , and compute ROM . 1

## 3: repeat

- 4: Phase 1: Predictor-Corrector
- 5: Set p i = p i -1 + u 0 . If p i &gt; p U , set p i = p U .
- 6: Build ROM cand i with a MOR algorithm.
- 7: if i=2 then

8:

$$\underset { \text{$\mathcal{P}_{\mathcal{OR}$} _ { i } } { \overline { R O M } } _ { i } = \mathcal { P } _ { \overline { R O M } _ { i - 1 } } ( \text{ROM} _ { i } ^ { \text{cand} } ).$$

9:

10:

11:

12:

13:

14:

15:

16:

17:

18:

19:

20:

21:

22:

23:

24:

25:

26:

27:

28:

29:

30:

31:

32:

## else

Compute ROM pred

i

of ROM

j

(

j

= 1 2

using extrapolation or regression

,

if

r

(

ROM

i

-

, . . . , i

1

,

-

1

).

ROM cand

i

)

&gt;

r

then

$$\underset { \alpha \colon } { \overline { R O M } } _ { i } = \mathcal { P } _ { R O M _ { i } ^ { \text{pred} } } ( \text{ROM} _ { i } ^ { \text{cand} } ).$$

else

$$\underset { \dots } { \overline { R O M } } _ { i } = \mathcal { P } _ { \overline { R O M } _ { i - 1 } } ( \text{ROM} _ { i } ^ { \text{cand} } ).$$

end if end if

## Phase 2: Adaptive Refinement

Set i h = i -1 .

## repeat

At p t = p i h + p i h +1 2 , build an interpolated ROM, namely ROM I,t , and a true ROM, namely ROM T,t . if e ( ROM I,t , ROM T,t ) &lt; τ e then

$$i _ { h } = i _ { h } + 1.$$

## else

$$\text{else} \quad \text{$i=i+1$.} \quad \text{$forj=i-1$ to $i_{h}+1$ do} \quad \overline { \text{ROM}_{j+1} \leftarrow \overline { \text{ROM}_{j}$, $p_{j+1} \leftarrow p_{j}$.} \quad \text{end for} \quad \overline { \text{ROM}_{i_{h}+1}=\mathcal{P}_{\overline{ROM}_{i_{h}}(\text{ROM}^{T,t)}$.} \quad \text{end for} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{endif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end of if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{start of the last $1$} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{endelse} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end endif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end endif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{endif} \quad \text{endif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{endif} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end elif} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{end else} \quad \text{start of the last $1$} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \ \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \ \text{end if} \quad \text{end ifif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{start if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{beginif} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{begin if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{ends if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{end if} \quad \text{in} \quad \text{end if} \quad \text{end if} \quad \text{end if}$$

end for end if

(

ROM pred

i

,

ROM cand

i

)

- 33: until i h = i .
- 34: until p i = p U .

## 35: Optional Post-Processing: Regression

- 36: for u = 1 to n d , v = 1 to 4 do
- 37: Apply the regression algorithm to D 1 ( u, v ) , D 2 ( u, v ) , . . . , D u,v i ( ) to compute the coefficients of the polynomial d r 0 ( u, v ) , d r 1 ( u, v ) , . . . , d r q ( u, v ) .

## 38: end for

- 39: for u = 1 to n s , v = 1 to 2 do
- 40: Apply the regression algorithm to S 1 ( u, v ) , S 2 ( u, v ) , . . . , S i ( u, v ) to compute the coefficients of the polynomial s r 0 ( u, v ) , s r 1 ( u, v ) , . . . , s r q ( u, v ) .

## 41: end for

Note that when pole-matching becomes difficult, i.e., too many swaps take place, we can simply give it up, save the ROM (with pole unmatched) for the moment and insert another ROM for refinement.The pole-matching will be done automatically at a later stage. A further remark is that in line 29, we can also conduct fewer pole-matchings, e.g., only for i h +1 .

## D. The Offline-Online Strategy

The proposed PMOR method can be implemented in an offline-online manner.

- · In the offline phase, we use Algorithm IV.1 to explore the parameter range of interest. It is relatively expensive: we need to adaptively choose the parameter values, at which we build ROMs, compute the ROMs, conduct pole matching, and optionally, perform regression.
- · In the online phase, we simply interpolate the repository ROMs or evaluate the polynomials computed from regression. The online phase is computationally very cheap.

## V. NUMERICAL RESULTS

In this section, we use two numerical examples to test the adaptive pole-matching PMOR method.

## A. The Nonlinear Parametric ' FOM ' Model

This academic example is adapted from the parametric ' FOM ' model in [19] to introduce more complex system dynamics such as nonlinear parametric dependency and pole crossing, since in the original model, the parameter dependency is linear . The model is of the form 3

$$( s \mathcal { I } - \mathcal { A } ( p ) ) X ( s, U ) & = \mathcal { B } U, \\ Y & = \mathcal { C } X ( s, U ),$$

where C = [100 100 100 100 100 100 100 100 1 , , , , , , , , , . . . , 1] ∈ R 1 × 1008 , B = C T , and A = diag {A 1 , A 2 , A 3 , A 4 , 1 2 3 , , , . . . , 1000 } ∈ R 1008 × 1008 with

$$\cdots, \, \dots, \, \dots \,, \, \dots \, \\ \mathcal { A } _ { 1 } & = \left [ \begin{array} { c c c } 4 p - 4 2 & 8 p + 2 0 0 \\ - 8 p - 2 0 0 & 4 p - 4 2 \end{array} \right ], \\ \mathcal { A } _ { 2 } & = \left [ \begin{array} { c c c } 2 p - 5 0 & p ^ { 2 } + 4 p + 2 1 0 \\ - p ^ { 2 } - 4 p - 2 1 0 & 2 p - 5 0 \end{array} \right ],$$

3 The parametric ' FOM ' model in [19] is adapted from the nonparametric ' FOM ' model in [20]

$$\mathcal { A } _ { 3 } & = \left [ \begin{array} { c c c } - 2 5 + p & 1 0 0 + p ^ { 2 } \\ - 1 0 0 - p ^ { 2 } & - 2 5 + p \end{array} \right ], \\ \mathcal { A } _ { 4 } & = \left [ \begin{array} { c c c } - 2 5 + 2 p & 1 5 0 - p ^ { 2 } \\ - 1 5 0 + p ^ { 2 } & - 2 5 + 2 p \end{array} \right ]. \\ \dots & \cdots \, \ell \quad \dots \, \cdots \, \cdots \,.$$

In this model of order 1008, the parameter range of interest is [ -10 10] , . In Algorithm IV.1, we set the initial conditions u 0 = π 3 , τ e = 0 001 . and q = 5 . For all figures, the relative error is computed as | ∫ 1000 1 H ( ω,p ) -H ω,p ( ) d ω | | ∫ 1000 1 H ( ω,p )d ω | .

(a) The Poles of Repository ROMs

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-f246f721-393f-49ad-8c79-68f8cdfdbc7e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-e1f142cd-851d-42f5-8a8e-2ad213a3c356.jpg)

Fig. 1. Interpolation Results

(a) The Relative Error

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-dcbf1c83-d91a-4379-aba0-7b1e490a0694.jpg)

Fig. 2. Regression Results

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-307c5238-4d33-4da1-8a22-935af92b53b2.jpg)

(b) The Relative Error

(b) The Regressive Polynomials

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-7ab081f6-8cf2-4889-ab66-7206761c0200.jpg)

(a) Interpolation Error

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-55b73dfa-f3c6-42ff-ab59-f66ca11c6f36.jpg)

(b) Regression Error

(c) The Regressive Polynomials

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-9d17544f-e2cb-4e03-9ed2-fffda6265894.jpg)

Fig. 3. Results without Adaptive Refinement

In Fig. 1, we show the results of interpolation, i.e., when we skip the post-processing for regression and use the repository

ROMs directly for interpolation. In total, 24 repository ROMs are built. In Fig. 2, we show the numerical results when the post-processing for regression is also performed. We observe that using order-5 polynomials for regression results in pretty accurate results. The data storage is reduced by a factor of 24 5 . In Fig. 3, we plot the results when we skip the adaptive refinement. In this case, 21 repository ROMs are built and we observe that some poles are wrongly matched and the regression follows wrong branches. This shows the importance of the adaptive refinement in guaranteeing accuracy, especially when the step length is large. The much better numerical result using adaptive refinement is only at the cost of 3 more repository ROMs.

## B. The Branchline Coupler: a Microwave Device

In this section, we apply our method to a branchline coupler model, which is a discretization of a time-harmonic Maxwell's equation. The full-order model is of the form

$$( \mathcal { K } ( \mu ) - s \mathcal { M } ) & X = B, \\ Y & = C X, \text{ \quad \quad } ( 3 1 )$$

where K ( µ ) = 1 µ K 0 , K 0 , M∈ R 33051 × 33051 , B ∈ R 33051 × 1 , C ∈ R 1 × 33051 and s = ω 2 . For detailed description of the model, we refer to [21], [22]. The repository ROMs of order 10 are built by a Krylov method [23], in which we use the interpolating expansion points s = 1 25 50 75 100 , , , , and order 2 for each of these interpolating expansion points. Since the reduced matrix of K ( µ i ) is nonsingular, the repository ROMs can be easily written into the form (1). To build a pROM for the parameter range µ ∈ [0 5 . , 0 9] . , Algorithm IV.1 built 18 repository ROMs under the tolerance τ e = 0 00001 . . Figure 4 shows that both the interpolated pROM and the regressed pROM capture the dynamics of the system well.

## VI. FURTHER DISCUSSIONS

Here are some further discussions.

- · The regression is successful only when all matchings are correct for dominant poles. The influence of mismatched poles on the interpolated pROM is only local, which is normally not too bad because the wrongly matched poles normally have similar positions and residues. For regression, however, the influence will be global because wrong branches are chosen. In this case, we actually try to regress a non-smooth function with polynomials, which is difficult. This phenomenon is shown in Fig. 3.
- · We would also like to emphasize that the interpolation approach is more general than the regression approach. The interpolation approach works well as long as the parameter dependence can be locally captured by the basis functions, e.g., using polynomials for exponential functions or trigonometric functions. The regression method, on the other hand, works well only for the cases that the parameter dependence can be globally captured by the basis functions. The interpolation approach even has potential to deal with eigenvalue bifurcation: even without any further considerations, the interpolation-based pROM

(b) The Relative Error of the Interpolated pROM

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-95c73a3a-c145-4e66-9f58-6addd7202fc0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-330f344f-059c-4e4b-8cb1-b15e3c59fd67.jpg)

(c) The Relative Error of the Regressed pROM (order 7 polynomials are used for regression)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00820v1-figure-232e8f8a-e37e-4ff0-9b47-595ea499183d.jpg)

Fig. 4. Numerical Results for the Branchline Example. Here, Y and y stand for the outputs of the FOM and the pROM, respectively.

is accurate as long as the parameter value is not too close to the bifurcation point. More detailed treatments of eigenvalue bifurcation will be future work. Therefore, we use regression as an optional post-processing technique. When the error between the interpolation-based pROM and the regression-based pROM is not sufficiently small, we discard the regression-based pROM and just use the interpolation-based pROM.

- · The proposed method only works when Assumption IV.2 holds. When A ( ) i in (1) is non-normal, its poles can be highly sensitive to perturbations of the parameters. We observed this phenomenon when we build order-22 Krylov-type ROMs for the branchline coupler model.
- · The proposed method can be easily extended to multipleinput multiple-output systems using the results in [9].

- · Stability can easily be preserved since our method deals directly with the poles. It is apparent that using linear interpolation, the stability will automatically be preserved since the linear interpolation of two poles in the left half-plane still lies in the left half-plane. Using other interpolation methods, although the stability is generally not automatically preserved, we can check whether the maximal real part of the interpolated (regressed) pole location is negative. If it is not, we can simply try another interpolation (regression) method and check it again. In the worst case, we can resort to linear interpolation for the parameter interval on which stability is difficult to be preserved with other interpolation/regression methods.

## VII. CONCLUSIONS

An adaptive parametric reduced-order modeling method based on interpolation/regression of adaptively built ROMs is proposed. For correct interpolation/regression, we convert all ROMs to the pole-residue realization and propose to use a branch and bound method to conduct pole-matching efficiently. To explore the parameter range efficiently and accurately, we propose to use a predictor-corrector technique and an adaptive refinement strategy. A regression method is proposed to further reduce data storage. The numerical results verify the high accuracy and robustness of the proposed method.

## REFERENCES

- [1] P. Benner, S. Gugercin, and K. Willcox, 'A survey of model reduction methods for parametric systems,' SIAM Review , vol. 57, no. 4, pp. 483531, 2015.
- [2] L. Feng, Y. Yue, N. Banagaaya, P. Meuris, W. Schoenmaker, and P. Benner, 'Parametric modeling and model order reduction for (electro-)thermal analysis of nanoelectronic structures,' J. Math. Ind. , vol. 6, no. 1, pp. 1-10, 2016.
- [3] Y. Yue and K. Meerbergen, 'Using Krylov-Pad´ e model order reduction for accelerating design optimization of structures and vibrations in the frequency domain,' International Journal for Numerical Methods in Engineering , vol. 90, no. 10, pp. 1207-1232, 2012.
- [4] M. Saadvandi, K. Meerbergen, and W. Desmet, 'Parametric dominant pole algorithm for parametric model order reduction,' Journal of Computational and Applied Mathematics , vol. 259, pp. 259-280, 2014.
- [5] Y. Zhang, L. Feng, S. Li, and P. Benner, 'Accelerating PDE constrained optimization by the reduced basis method: application to batch chromatography,' Internat. J. Numer. Methods Engrg. , vol. 104, no. 11, pp. 983-1007, 2015.
- [6] G. Rozza, D. B. P. Huynh, and A. T. Patera, 'Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations,' Archives of Computational Methods in Engineering , vol. 15, no. 3, pp. 229-275, 2008.
- [7] U. Baur, C. A. Beattie, P. Benner, and S. Gugercin, 'Interpolatory projection methods for parameterized model reduction,' SIAM J. Sci. Comput. , vol. 33, no. 5, pp. 2489-2518, 2011.
- [8] D. Amsallem and C. Farhat, 'Interpolation method for the adaptation of reduced-order models to parameter changes and its application to aeroelasticity,' AIAA J. , vol. 46, pp. 1803-1813, 2008.
- [9] Y. Yue, L. Feng, and P. Benner, 'Reduced-order modelling of parametric systems via interpolation of heterogeneous surrogates,' Advanced Modeling and Simulation in Engineering Sciences , pp. 1-28, 2019, to appear.
- [10] H. Panzer, J. Mohring, R. Eid, and B. Lohmann, 'Parametric model order reduction by matrix interpolation,' at-Automatisierungstechnik , vol. 58, no. 8, pp. 475-484, 2010.
- [11] D. Amsallem and C. Farhat, 'An online method for interpolating linear parametric reduced-order models,' SIAM J. Sci. Comput. , vol. 33, no. 5, pp. 2169-2198, 2011.
- [12] U. Baur and P. Benner, 'Modellreduktion f¨ ur parametrisierte Systeme durch balanciertes Abschneiden und Interpolation (Model Reduction for Parametric Systems Using Balanced Truncation and Interpolation),' atAutomatisierungstechnik , vol. 57, no. 8, pp. 411-420, 2009.
- [13] U. Baur, P. Benner, A. Greiner, J. G. Korvink, J. Lienemann, and C. Moosmann, 'Parameter preserving model reduction for MEMS applications,' Math. Comput. Model. Dyn. Syst. , vol. 17, no. 4, pp. 297317, 2011.
- [14] Y. Yue, L. Feng, and P. Benner, 'Interpolation of reduced-order models based on modal analysis,' in 2018 IEEE MTT-S International Conference on Numerical Electromagnetic and Metaphysics Modeling and Optimization (NEMO) , 2018.
- [15] --,'An adaptive method for interpolating reduced-order models based on matching and continuation of poles,' in 2019 IEEE MTT-S International Conference on Numerical Electromagnetic and Metaphysics Modeling and Optimization (NEMO) , 2019.
- [16] G. H. Golub and C. F. van Van Loan, Matrix Computations , 3rd ed. Baltimore and London: The Johns Hopkins University Press, 1996.
- [17] J. Clausen, 'Branch and bound algorithms - principles and examples,' University of Copenhagen, Tech. Rep., 1999.
- [18] N. Martins, L. T. G. Lima, and H. J. C. P. Pinto, 'Computing dominant poles of power system transfer functions,' IEEE Trans. Power Syst. , vol. 11, no. 1, pp. 162-170, 1996.
- [19] A. C. Ionita and A. C. Antoulas, 'Data-driven parametrized model reduction in the Loewner framework,' SIAM J. Sci. Comput. , vol. 36, no. 3, pp. A984-A1007, 2014.
- [20] Y. Chahlaoui and P. Van Dooren, 'A collection of benchmark examples for model reduction of linear time invariant dynamical systems,' SLICOT Working Note, Tech. Rep. 2002-2, 2002, available from www.slicot.org.
- [21] The MORwiki Community, 'Branchline coupler,' MORwiki -Model Order Reduction Wiki, 2013. [Online]. Available: http: //modelreduction.org/index.php/Branchline Coupler
- [22] M. W. Hess and P. Benner, 'Fast evaluation of time-harmonic Maxwell's equations using the reduced basis method,' IEEE Trans. Microw. Theory Techn. , vol. 61, no. 6, pp. 2265-2274, 2013.
- [23] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems , ser. Advances in Design and Control. Philadelphia, PA: SIAM Publications, 2005, vol. 6.
Yao Yue Lihong Feng Peter Benner
An Adaptive Pole-Matching Method for Interpolating Reduced-Order Models
An adaptive parametric reduced-order modeling method based on interpolating
poles of reduced-order models is proposed in this paper. To guarantee correct
interpolation, a pole-matching process is conducted to determine which poles of
two reduced-order models correspond to the same parametric pole. First, the
pole-matching in the scenario of parameter perturbation is discussed. It is
formulated as a combinatorial optimization problem and solved by a branch and
bound algorithm. Then, an adaptive framework is proposed to build repository
ROMs at adaptively chosen parameter values, which well represent the parameter
domain of interest. To achieve this, we propose techniques including a
predictor-corrector strategy and an adaptive refinement strategy, which enable
us to use larger steps to explore the parameter domain of interest with good
accuracy. The framework also consists of regression as an optional
post-processing phase to further reduce the data storage. The advantages over
other parametric reduced-order modeling approaches are, e.g., compatibility
with any model order reduction method, constant size of the parametric
reduced-order model with respect to the number of parameters, and capability to
deal with complicated parameter dependency.