# Comparison of parametric model order reduction methods to solve magneto-quasistatic and electro-quasistatic problems

Wei Chen, Thomas Henneron, Stephane Clenet

► To cite this version:

Wei Chen, Thomas Henneron, Stephane Clenet. Comparison of parametric model order reduction methods to solve magneto-quasistatic and electro-quasistatic problems. Finite Elements in Analysis and Design, 2025, 252, 10.1016/j.finel.2025.104444. hal-05262715

HAL Id: hal-05262715

https://hal.science/hal-05262715v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

# Comparison of parametric model order reduction methods to solve magneto-quasistatic and electro-quasistatic problems

Wei Chen, Thomas Henneron, Stéphane Clénet*

Univ. Lille, Arts et Metiers Institute of Technology, Centrale Lille, Junia, ULR 2697 - L2EP, Lille, 59000, France

## ARTICLE INFO

**Keywords:**

Model Order Reduction

Finite Element method

Electro-quasistatic problem

Magneto-quasistatic problem

Parameterized Model Order Reduction

Multi-moment matching method

Interpolation of projection subspaces method

## ABSTRACT

In this paper, we compare two parametric model order reduction methods, the multi-moment matching method and the interpolation of projection subspaces method for the magneto-quasistatic (MQS) and electro-quasistatic (EQS) problems derived from Maxwell's equations and discretized with the Finite Element (FE) method. The two problems considered are both governed by the differential-algebraic equations. The material characteristic parameters as well as the geometry parameters have been considered. The applications are two realistic test cases: an EQS model of a transformer bushing under insulation defect uncertainty and a MQS model of a planar inductor with geometric and material variations. The result shows that both methods approximate well global quantities, such as the current or the voltage, as well as the local quantities like field distributions. The multi-moment matching method remains always faster in the online stage, since the reduced basis is not parameter dependent, requiring no reduced basis calculation. The multi-moment matching method requires an affine decomposition of the FE model, which is not easy to obtain when considering geometry parameters. A hybrid method is proposed and tested leading to more accurate results than the interpolation of projection subspaces method but much easier to implement than the multi-moment matching method.

## 1. Introduction

For the optimization and fault detection during and after the design stages of electromagnetic devices, numerical models based on the Finite Element (FE) method are usually considered as the reference model due to their high fidelity in the field of low-frequency computational electromagnetism. These numerical models are based on Maxwell's equations. When operating at low frequencies, Maxwell's equations can be simplified into two quasistatic problems: the magneto-quasistatic (MQS) problem and the electro-quasistatic (EQS) problem. These simplifications make the problems easier to solve with the FE method without significantly losing precision. Applying the FE method to solve the MQS and EQS problems often results in a set of Differential-Algebraic Equations (DAE).

In practice, the FE model has a large number of unknowns due to the need to discretize complex geometric shapes accurately. This leads to significant computation times, which prevent the intensive use of such models. Different Model Order Reduction (MOR) methods have been proposed to avoid repeatedly solving these large-size FE models, also known as full models. Rather than solving the full model, MOR methods aim to build a reduced model through projection of the full one onto an appropriate basis, which can be constructed using methods such as the Proper Orthogonal Decomposition (POD) method [1] or the Krylov subspace method [2] for example.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: stephane.clenet@ensam.eu (S. Clénet).</footnote>
Often, the dimensions of the geometry, the material characteristics, or the sources are not perfectly known and are represented by parameters. In that case, the FE model will lead to a parameterized equation system. This system must be solved for each new set of parameters, leading to a significant computational time. MOR methods can be applied to construct a reduced basis for a given set of parameters using the POD method or the Krylov subspace approach, for example. This leads to a reduced model that is generally not valid for the entire parameter space. It is necessary to adapt the reduced basis according to the set of parameters considered. A parametric order reduction method (PMOR) should be used [3–7]. The PMOR method can result in a reduced basis or a reduced model which remains a function of the set of parameters. The PMOR method greatly speeds up the construction of the reduced basis when the parameters change. The PMOR method can be used advantageously in various applications, such as design optimization or uncertainty quantification, for which the full model needs to be solved multiple times for different sets of parameters.

In this paper, we are interested in the application of PMOR methods to reduce MQS and EQS problems discretized by the FE method and parameterized by d parameters $\mathbf{p} = [p_1, ..., p_d] \in \mathbb{R}^d$. These problems can be written under the general form:

$$ 
\begin{aligned}
\mathbf{K}(\mathbf{p})\mathbf{X}(s, \mathbf{p}) + s\mathbf{N}(\mathbf{p})\mathbf{X}(s, \mathbf{p}) &= U(s)\mathbf{F}(\mathbf{p}), \\
y(s, \mathbf{p}) &= \mathbf{F}(\mathbf{p})^T\mathbf{X}(s, \mathbf{p}).
\end{aligned} 
\quad (1) $$

Where $\mathbf{K}(\mathbf{p})$, $\mathbf{N}(\mathbf{p}) \in \mathbb{R}^{N \times N}$, $\mathbf{F}(\mathbf{p}) \in \mathbb{R}^{N \times 1}$ are state-space matrices and s is the Laplace variable. We denote by $\mathbf{X}(s, \mathbf{p}) \in \mathbb{R}^{N \times 1}$ the state vector and $y(s, \mathbf{p})$ the quantity of interest.

The PMOR method aims to find an approximation of $\mathbf{X}(s, \mathbf{p})$ into a subspace $\mathcal{V}(\mathbf{p})$, spanned by the reduced basis $\mathbf{V}(\mathbf{p})$, and project the system on the test subspace $\mathcal{W}(\mathbf{p})$, spanned by the basis $\mathbf{W}(\mathbf{p})$. The dimension r of $\mathcal{V}(\mathbf{p})$ is much smaller than the size of the full model (1) ($r \ll N$). In this paper, only the Ritz-Galerkin projection method will be considered that is to say, the reduced bases are the same $\mathbf{V}(\mathbf{p}) = \mathbf{W}(\mathbf{p})$. PMOR methods aim then to construct the reduced basis $\mathbf{V}(\mathbf{p})$ such that the solution $\mathbf{X}(s, \mathbf{p})$ can be approximated by:

$$ \mathbf{X}(s, \mathbf{p}) \approx \mathbf{V}(\mathbf{p})\mathbf{X}_r(s, \mathbf{p}) \quad (2) $$

where $\mathbf{X}_r(s, \mathbf{p}) \in \mathbb{R}^{N \times 1}$ is the reduced state vector whose dimension r is much lower than the dimension N of $\mathbf{X}(s, \mathbf{p})$. It should be mentioned that, depending on the PMOR method, the reduced basis $\mathbf{V}(\mathbf{p})$ can depend or not on the parameter vector $\mathbf{p}$. The reduced state vector satisfies the following equation system after applying a Galerkin projection method:

$$ 
\begin{aligned}
\mathbf{V}(\mathbf{p})^T \mathbf{K}(\mathbf{p}) \mathbf{V}(\mathbf{p}) \mathbf{X}_r(s, \mathbf{p}) + s \mathbf{V}(\mathbf{p})^T \mathbf{N}(\mathbf{p}) \mathbf{V}(\mathbf{p}) \mathbf{X}_r(s, \mathbf{p}) &= U(s) \mathbf{V}(\mathbf{p})^T \mathbf{F}(\mathbf{p}), \\
y(s, \mathbf{p}) &= \mathbf{F}^T(\mathbf{p}) \mathbf{V}(\mathbf{p}) \mathbf{X}_r(s, \mathbf{p}).
\end{aligned} 
\quad (3) $$

Which leads to the reduced system of equations:

$$ 
\begin{aligned}
\mathbf{K}_r(\mathbf{p})\mathbf{X}_r(s, \mathbf{p}) + s\mathbf{N}_r(\mathbf{p})\mathbf{X}_r(s, \mathbf{p}) &= U(s)\mathbf{F}_r(\mathbf{p}), \\
y(s, \mathbf{p}) &= \mathbf{F}_r^T(\mathbf{p})\mathbf{X}_r(s, \mathbf{p}).
\end{aligned} 
\quad (4) $$

With the reduced matrices and vector given by:

$$ \mathbf{K}_r(\mathbf{p}) = \mathbf{V}(\mathbf{p})^T \mathbf{K}(\mathbf{p}) \mathbf{V}(\mathbf{p}), \quad \mathbf{N}_r(\mathbf{p}) = \mathbf{V}(\mathbf{p})^T \mathbf{N}(\mathbf{p}) \mathbf{V}(\mathbf{p}), \quad \mathbf{F}_r(\mathbf{p}) = \mathbf{V}(\mathbf{p})^T \mathbf{F}(\mathbf{p}). \quad (5) $$

The size $r \times r$ of the reduced system of Eqs. (4) is much smaller than the size $N \times N$ of the initial problem (1) which leads to a gain in terms of computational time and also in terms of memory usage, as only the solutions of the reduced system need to be stored rather than those of the initial problem.

The projection-based PMOR methods, applied to construct the reduced basis $\mathbf{V}(\mathbf{p})$, can be classified into four categories. The first category of methods aims at constructing one global basis $\mathbf{V}$ valid for the entire parameter space. This global basis could be formed by the concatenation of several bases, each spanning some Krylov subspaces. These Krylov subspaces [8–10] implicitly match some generalized moments of the transfer function $H(s, \mathbf{p}) = y(s, \mathbf{p})/U(s)$. Those methods are the so-called Multi-(Parameter) Moment Matching (MultiPMomMtch) or multivariate Padé approximation [4]. The balanced-truncation method [11] can also be used to construct such a global basis.

The second category of methods is based on the construction of several local bases $\mathbf{V}(\mathbf{p})$ at several sample points of the parameter space $\mathbf{p}_1, ..., \mathbf{p}_K$. We denote by $\mathbf{V}_k$ the corresponding local bases obtained for the parameter $\mathbf{p}_k$, $k = 1, ..., K$. For a new parameter $\mathbf{p}$, a reduced base $\mathbf{V}(\mathbf{p})$ is obtained by interpolating the subspaces spanned by $\mathbf{V}_k$, with $k = 1, ..., K$, in a tangent space to a manifold of these subspaces [12]. This method is so-called interpolation of projection subspaces.

The third category of methods interpolates directly the local reduced-order coefficient matrices $\mathbf{K}_{r,k}(\mathbf{p}) = \mathbf{V}_k^T \mathbf{K}(\mathbf{p}) \mathbf{V}_k$ and $\mathbf{N}_{r,k}(\mathbf{p}) = \mathbf{V}_k^T \mathbf{N}(\mathbf{p}) \mathbf{V}_k$ corresponding to the parameters $\mathbf{p}_1, ..., \mathbf{p}_K$ by using matrix manifold interpolation [13]. This method is so-called interpolation in the time domain [14].

The fourth category of methods interpolates the local transfer functions $H_{r,k}(s, \mathbf{p})$ constructed using the local basis $\mathbf{V}_k$ [15,16]. Some different PMOR methods, except for the second category of methods, have been compared on examples [4] based on systems of ordinary differential equations (ODE).

The fourth category of methods focuses on the calculation of the quantity of interest $y(s, \mathbf{p})$, which is often a global quantity, but does not lead to any information about local quantities, like field distribution, reconstructed from the true or approximated $\mathbf{X}$. In practice, in computational electromagnetics, we are interested in local quantities of interest, such that field distributions
for example. Indeed, these distributions enable to have access to the distribution of losses, enabling to evaluate the temperature distribution, quantities important to control in order to avoid any failure or quick aging of an electromagnetic device.

The third category of methods, interpolation in the time domain, is proposed for ordinary differential equations. For MQS and EQS problems, the system of equations is rather differential-algebraic equations (DAE). Indeed, the matrix N is singular. Hence, the third category of methods cannot be applied straightforwardly since an appropriate manifold related to the reduced matrix N, cannot be found.

To reduce our parametric EQS or MQS problems, we can only apply directly the PMOR methods of the first and second categories, so in the following, we will focus mainly on the application of these methods to reduce EQS and MQS problems.

The MultiPMomMtch method has been applied to electromagnetic problems discretized by the integral equation method [17], by the Partial Element Equivalent Circuit (PEEC) method [18,19] and by the Finite Element method in computational electromagnetics at high frequency (full wave) [20,21]. Specifically, in [20], the MultiPMomMtch method was first used to construct a parametric reduced model considering two material characteristics and the frequency as parameters. The MultiPMomMtch method has also been used for geometric parametrization in mechanics [22] using the isoparametric Finite Element method. The interpolation of projection subspaces method has been applied to solve the MQS problem using the Finite Element method [23] for the non linear problem. However, the interpolation of projection subspaces and the MultiPMomMtch methods have not been used to study MQS and EQS problems with geometric or material characteristic parameters.

Moreover, as mentioned before, the comparison of the MultiPMomMtch method and the interpolation of projection subspaces method has not been done in the context of the DAE problem and more especially in the context of low frequency computational electromagnetism. Finally, it seems that the hybrid method combining both methods has not been studied yet.

In this paper, we will apply the two PMOR methods, MultiPMomMtch and the interpolation of projection subspaces methods, to reduce parametric FE EQS and MQS problems. Parameters related to geometry and material characteristics will be considered simultaneously. The two methods will be compared and a hybrid method combining both will be proposed.

The paper is organized as follow, in Section 2, the MultiPMomMtch method and the interpolation of projection subspaces method will be presented. In Section 3, the parametric FE EQS and MQS problems will be described. Then the two PMOR methods will be applied to the EQS and MQS problems. For the EQS problem, a resin-impregnated paper bushing of a high-voltage transformer will be considered. The effect of a default of insulation will be analyzed. In that case, the problem becomes parametric since the permittivity and conductivity of the defect are unknowns. For the MQS problem, we will consider a planar inductor used in power electronic converters. In that case, we are interested in the influence of the size of the airgap and the position of the conductor as well as the permeability of the ferrite core. In both cases, the two PMOR methods will be compared in terms of accuracy, computational time and implementation. Finally, a hybrid method combining both methods will be developed and applied in the case of the planar inductor. This hybrid method will also be compared with the two original methods.

## 2. PMOR methods

### 2.1. Multi-moment matching method

At the origin, the multi-moment matching method [9,24] was proposed to study large interconnected networks. The multi-moment matching methods were introduced to take into account the effect of parameters related to material characteristics, sources and geometry. Similarly to the Passive Reduced-order Interconnect Macromodeling Algorithm (PRIMA) method [2], it is shown [8] that the explicit calculation of moments often results in poor precision compared to the implicit calculation of the moment, which uses, as an intermediate, a reduced basis that needs to be determined. It should be mentioned that this method will lead to a reduced basis V that does not depend on the parameters p.

In the following, we will present the algorithm to determine the reduced basis V, the global basis, [8] for the single input single output (SISO) case, providing that it can be easily generalized to the multi input multi output (MIMO) case. Considering a parametric equation system (1), the aim of the PMOR via the implicit moment matching method is to construct a reduced basis V and use the Galerkin projection method, which transforms the full model (1) into (3).

The transfer function $H(s, p) = y(s, p)/U(s)$ defined from (1) can be rewritten under the form:

$$H(s, p) = H(\hat{p}) = F^T(p)Q^{-1}(\hat{p})F(p) \qquad (6)$$

where $Q(\hat{p}) = K(p) + sN(p)$ and $\hat{p} := (s, p) \in \mathbb{C}^{1 \times (d+1)}$.

To apply the MultiPMomMtch method, it is required that the system matrix $Q(\hat{p})$ can be written in the following affine form:

$$Q(\hat{p}) = Q_0 + h_1(\hat{p})Q_1 + \dots + h_p(\hat{p})Q_p \qquad (7)$$

where $h_i(\hat{p}) : \mathbb{C}^{d+1} \to \mathbb{C}$ are the scalar functions of $\hat{p}$ and $p$ is the number of terms of the decomposition (7).

The Taylor series expansion of the state variable $X(\hat{p})$ at an expansion point $\hat{p}_i := (s_i, p_i)$ is:

$$
\begin{aligned}
X(\hat{p}) &= [Q(\hat{p})]^{-1}F(p)U(s) \\
&= [Q_0 + h_1(\hat{p})Q_1 + \dots + h_p(\hat{p})Q_p]^{-1}F(p)U(s) \\
&= [Q(\hat{p}_i) + \sigma_1 Q_1 + \dots + \sigma_p Q_p]^{-1}F(p)U(s) \\
&= [I - (\sigma_1 M_1 + \dots + \sigma_p M_p)]^{-1}F_M(p)U(s) \\
&= \sum_{i=0}^{\infty} [\sigma_1 M_1 + \dots + \sigma_p M_p]' F_M(p)U(s)
\end{aligned}
\qquad (8)
$$
<table><caption>**Table 1**
Synthesis of parameters for the MultiPMomMtch method.</caption><thead><tr><th>d</th><th>p</th><th>q</th><th>l</th></tr></thead><tbody><tr><td>Number of parameters</td><td>Number of affine functions</td><td>Order of series expansion</td><td>Number of expansion points</td></tr></tbody></table>

where $\sigma_j = h_j(\hat{p}) - h_j(\hat{p}_j)$ and $M_j = -[Q(\hat{p}_j)]^{-1}Q_j$, $j = 1, 2, \dots, p$. In the SISO case, the vector $F_M(p)$ is defined as:

$$F_M(p) = [Q(\hat{p}_i)]^{-1}F(p) \qquad (9)$$

Remark: the matrix $Q(\hat{p}_i)$ is in general not invertible for the MQS problem. However, if the source term $F_x$ is, the equation system $Q(\hat{p}_i)X = F_x$ can be solved by an iterative solver, the system being then implicitly autogauged [25].

In the general case where $F(p)$ depends on the parameter, $F(p)$ needs to be approximated by an affine form similar to (7). In our applications, the vector $F(p)$ is independent of the parameters $p$, that is to say, $F(\hat{p}) = F$. The vector $F_M$ is independent of parameter $p$ once the expansion point $\hat{p}$ is fixed.

By defining

$$ \begin{aligned} x_0 &= F_M, \\ x_1 &= [\sigma_1 M_1 + \dots + \sigma_p M_p] F_M, \\ x_2 &= [\sigma_1 M_1 + \dots + \sigma_p M_p]^2 F_M, \\ \vdots \\ x_q &= [\sigma_1 M_1 + \dots + \sigma_p M_p]^q F_M, \\ \vdots \end{aligned} \qquad (10) $$

We have $X = (x_0 + x_1 + \dots + x_q + \dots)U(s)$ and obtain the following recursive relationship:

$$ \begin{aligned} x_0 &= F_M, \\ x_1 &= [\sigma_1 M_1 + \dots + \sigma_p M_p]x_0, \\ x_2 &= [\sigma_1 M_1 + \dots + \sigma_p M_p]x_1, \\ \vdots \\ x_q &= [\sigma_1 M_1 + \dots + \sigma_p M_p]x_{q-1}, \\ \vdots \end{aligned} \qquad (11) $$

We can then define a sequence of calculations of matrices $R_q$ belonging to $\mathbb{R}^{N \times p^q}$ such that:

$$ \begin{aligned} \mathbf{R}_0 &= \mathbf{F}_M, \\ \mathbf{R}_1 &= [\mathbf{M}_1 \mathbf{R}_0, \mathbf{M}_2 \mathbf{R}_0, \dots, \mathbf{M}_p \mathbf{R}_0], \\ \vdots \\ \mathbf{R}_q &= [\mathbf{M}_1 \mathbf{R}_{q-1}, \mathbf{M}_2 \mathbf{R}_{q-1}, \dots, \mathbf{M}_p \mathbf{R}_{q-1}], \\ \vdots \end{aligned} \qquad (12) $$

According to the expression of the state vector $X(\hat{p})$ in (8), it can be shown that this vector belongs to the space defined by $\text{span}\{\mathbf{R}_0, \dots, \mathbf{R}_j, \dots\}$. By truncating the space up to the $q$th first terms, we can obtain the reduced basis $\mathbf{V}_{\hat{p}_i}$ defined as:

$$ \text{range}(\mathbf{V}_{\hat{p}_i}) = \text{span}\{\mathbf{R}_0, \dots, \mathbf{R}_j, \dots, \mathbf{R}_q\} \qquad (13) $$

It should be noted that the number of columns in the matrix $\mathbf{R}_j \in \mathbb{R}^{N \times p^j}$ increases exponentially with $j$. It is usually recommended to set $q \le 1$ and use only the values $\mathbf{R}_0$ and $\mathbf{R}_1$ [26]. The matrix $\mathbf{V}_{\hat{p}_i}$ is constructed at the expansion point $\hat{p}_i$. The global basis matrix is obtained by concatenation and orthogonalization of $l$ matrices computed at $l$ expansion points: $\hat{p}_1, \hat{p}_2, \dots, \hat{p}_l$.

$$ \mathbf{V} = \text{orth}\{\mathbf{V}_{\hat{p}_1}, \mathbf{V}_{\hat{p}_2}, \dots, \mathbf{V}_{\hat{p}_l}\} \qquad (14) $$

The synthesis of different integer parameters for the MultiPMomMtch method is given in Table 1. For $i = 1, \dots, l$, the matrix $\mathbf{V}_{\hat{p}_i}$ is computed using Algorithm 1 by initializing the matrix $\hat{\mathbf{Q}}$ as:

$$ \hat{\mathbf{Q}} = \mathbf{Q}(\hat{\mathbf{p}}_i) \qquad (15) $$

Since the matrices $\mathbf{Q}_1, \dots, \mathbf{Q}_p$ are constant, it means that Algorithm 1 in Appendix A is called $l$ times changing only the matrix $\hat{\mathbf{Q}}$ defined in (15) as input. Algorithm 1 can be understood as a repeated modified Gram-Schmidt (MGS) process [8]. In Step 20, $\|\mathbf{w}\|_2 > \epsilon$ is the criterion used to avoid deflation of vectors in $\mathbf{R}_j$. If $\|\mathbf{w}\|_2 = 0$, it means that $\mathbf{w}$ is not linearly independent of the previously computed vectors $\mathbf{v}_1, \mathbf{v}_2, \dots, \mathbf{v}_{\text{col}-1}$ with $\text{col}$ defined in Step 16. In this case, the vector $\mathbf{w}$ is dropped (the detection of this case is done in Step 20).
## 2.2. Interpolation of projection subspaces

In this section, the method of interpolation of projection subspaces will be reviewed. We will first present some important properties of the Grassmann manifold. However, we will not recall all underlying theories of differential geometry. The reader can find an introduction in [27–29].

Let $r$ and $N$ be positive integers with $r \le N$. The Grassmann manifold is the set of $r$-dimensional subspaces of $\mathbb{R}^N$, denoted here $G_r(\mathbb{R}^N)$. More details about the Grassmann manifolds can be found in [27,29]. Let $\mathbf{Y}$ be a matrix of $\mathbb{R}^{r \times N}$ of rank $r$, We denote $y_1, \dots, y_r$ the columns of $\mathbf{Y}$. The linear subspace $\mathcal{Y}$ spanned by the $y_1, \dots, y_r$ belongs to $G_r(\mathbb{R}^N)$:

$$ \mathcal{Y} = \text{range}(\mathbf{Y}) \in G_r(\mathbb{R}^N) \qquad (16) $$

The matrix $\mathbf{Y}$ is called as a matricial representation of the subspace $\mathcal{Y}$. For $\mathbf{M}$ in general linear group of degree $r$ $GL(r)$, the set of $r \times r$ non singular matrices, the range of the matrix $\mathbf{YM}$ is the same as the one of $\mathbf{Y}$ that is to say $\mathcal{Y}$:

$$ \mathcal{Y} = \text{range}(\mathbf{YM}) = \text{range}(\mathbf{Y}) \qquad (17) $$

If we consider $K$ parameter values $\mathbf{p}_j$, $j = 1, \dots, K$, and the associated reduced basis $\mathbf{W}_j$ calculated using POD (This matrix is calculated from snapshots solving the equation system (1) for different values of $s$). We want now to determine the reduced basis $\mathbf{V}(\mathbf{p})$, the local basis, for any $\mathbf{p}$ of the parameter space. It is shown in [12] that a classical interpolation of $\mathbf{V}(\mathbf{p})$ from $\mathbf{p}_j$, $j = 1, \dots, K$, such as the interpolation of Lagrangian-type matrices, may not always result in a matrix of full rank. To solve this problem, it is then proposed to carry out the interpolation on the tangent space of a point of the Grassmann manifold. This process is composed of four steps [12,14], which will be presented in the following.

We denote $T_{W_0} G_r(\mathbb{R}^N)$ the tangent space of $G_r(\mathbb{R}^N)$ at $W_0 \in G_r(\mathbb{R}^N)$, $\text{Log}_{W_0}$ the logarithm of $W_0 \in G_r(\mathbb{R}^N)$ and $\text{Exp}_{W_0}(\mathcal{Y})$ the exponent of $\mathcal{Y} \in T_{W_0} G_r(\mathbb{R}^N)$. The exponential mapping $\text{log}_{W_0}(\mathbf{W})$ is defined via a geodesic path on $G_r(\mathbb{R}^N)$. A geodesic path connecting two points of a manifold is the “locally shortest” curve that relates $W_0$ and $\mathbf{W}$. For a point $W_0 \in G_r(\mathbb{R}^N)$, the logarithm mapping $\text{Log}_{W_0}$ is the inverse map of $\text{Exp}_{W_0}(\mathcal{Y})$ in the neighborhood of $W_0$. More details can be found in [27].

We will present the interpolation process from the $K$ subspaces $W_1 \in G_r(\mathbb{R}^N), \dots, W_K \in G_r(\mathbb{R}^N)$ with a matricial representation $\mathbf{W}_1, \dots, \mathbf{W}_K$. In our case, the matrices $\mathbf{W}_j$, $j = 1, \dots, K$, are the POD basis matrices, computed from the solutions of the FE model at different frequencies, for the parameter $\mathbf{p}_j$. The $K$ local basis matrices $\mathbf{W}_1, \dots, \mathbf{W}_K$ are precomputed in the offline stage. We now want to interpolate a new local basis matrix $\mathbf{V}(\mathbf{p})$ for a new set of parameters $\mathbf{p}$, which is spanned by the interpolated subspace $\mathcal{W}(\mathbf{p})$. The procedure is described below:

1. Choose a subspace $W_0$ among $W_1 \in G_r(\mathbb{R}^N), \dots, W_K \in G_r(\mathbb{R}^N)$ to construct the tangent space. This subspace is so-called contact point in the following.

2. Map $W_1 \in G_r(\mathbb{R}^N), \dots, W_K \in G_r(\mathbb{R}^N)$ to the tangent space $T_{W_0} G_r(\mathbb{R}^N)$ at the contact point by $\text{log}_{W_0}$. In order to do so, we need to compute the thin Singular Value Decompositions (SVD):

$$ (\mathbf{I} - \mathbf{W}_0(\mathbf{W}_0^T \mathbf{W}_0)^{-1} \mathbf{W}_0^T) \mathbf{W}_j (\mathbf{W}_0^T \mathbf{W}_j)^{-1} (\mathbf{W}_0^T \mathbf{W}_0)^{1/2} = \mathbf{U}_j \mathbf{\Sigma}_j \mathbf{V}_j^T \qquad (18) $$

where $j = 1, \dots, K$ and $j \neq i$. The image of $\mathcal{Y}_{W_j} = \text{Log}_{W_0}(W_j)$ is spanned by the columns of a matrix

$$ \mathbf{Y}_{W_j} = \mathbf{U}_j \arctan(\mathbf{\Sigma}_j) \mathbf{V}_j^T, \quad j = 1, \dots, K, j \neq i \qquad (19) $$

and $\mathcal{Y}_{W_0} = \text{Log}_{W_0}(W_0) = \mathbf{0}$.

3. Interpolate on $T_{W_0} G_r(\mathbb{R}^N)$ using a standard interpolation technique. For a parameter $\mathbf{p} \in \mathbb{R}^d$, compute the matrix

$$ \mathbf{Y}_{\mathbf{W}}(\mathbf{p}) = \sum_{j=1}^{K} f_j(\mathbf{p}) \mathbf{Y}_{W_j} \qquad (20) $$

where the functions $f_j(\mathbf{p})$ are the weight computed with the chosen interpolation method. The matrix $\mathbf{Y}_{\mathbf{W}}(\mathbf{p})$ spans the subspace $\mathcal{Y}_{\mathbf{W}}(\mathbf{p}) \in T_{W_0} G_r(\mathbb{R}^N)$. The weight function can be chosen under various forms. In this paper, multi linear interpolation functions are used. For example, let consider the case of dimension $d = 2$ and with four sample points lie on a grid: $\mathbf{p}_1 = (x_1, y_1)$, $\mathbf{p}_2 = (x_1, y_2)$, $\mathbf{p}_3 = (x_2, y_1)$ and $\mathbf{p}_4 = (x_2, y_2)$. The weight function at parameter $\mathbf{p} = (x, y)$ is then given by:

$$ f_1(\mathbf{p}) = \frac{x - x_2}{x_1 - x_2} \cdot \frac{y - y_2}{y_1 - y_2}, \quad f_2(\mathbf{p}) = \frac{x - x_2}{x_1 - x_2} \cdot \frac{y - y_1}{y_2 - y_1}, \qquad (21) $$

$$ f_3(\mathbf{p}) = \frac{x - x_1}{x_2 - x_1} \cdot \frac{y - y_2}{y_1 - y_2}, \quad f_4(\mathbf{p}) = \frac{x - x_1}{x_2 - x_1} \cdot \frac{y - y_1}{y_2 - y_1}. \qquad (22) $$

4. Map $\mathcal{Y}_{\mathbf{W}}(\mathbf{p})$ back to the Grassmann manifold $G_r(\mathbb{R}^N)$ by $\text{Log}_{W_0}$. In order to do so, we need first compute the thin SVD:

$$ \mathbf{Y}_{\mathbf{W}}(\mathbf{p}) = \mathbf{U}_{\mathbf{p}} \mathbf{\Sigma}_{\mathbf{p}} \mathbf{V}_{\mathbf{p}}^T \qquad (23) $$

and a matricial representation:

$$ \mathbf{W}(\mathbf{p}) = \mathbf{W}_0(\mathbf{W}_0^T \mathbf{W}_0)^{-1/2} \mathbf{V}_{\mathbf{p}} \cos(\mathbf{\Sigma}_{\mathbf{p}}) + \mathbf{U}_{\mathbf{p}} \sin(\mathbf{\Sigma}_{\mathbf{p}}) \qquad (24) $$

of the desired interpolated subspace $\mathcal{W}(\mathbf{p}) = \text{Exp}_{W_0}(\mathcal{Y}_{\mathbf{W}(\mathbf{p})})$.
By applying this algorithm, we can calculate a local basis matrix $V(p)$ for any parameter set $p$. If we divide the process into the offline and online stages, we have to calculate:

1. the reduced bases $V_k$ the corresponding local bases obtained for the parameter $p_k$, $k = 1, \dots, K$ from the full model (1) by applying the POD method in the offline stage,

2. an approximation of the reduced basis $V(p)$ by applying the interpolation described above in the online stage.

It exists also some acceleration techniques to make the interpolation of projection subspaces method real-time [14]. It is known that the interpolation of projection subspaces method is influenced by the chosen contact point $W_0$. There are some techniques that decrease the influence of the contact point. However, those techniques [30] would increase the number of SVD procedures, which would then increase the computational time of the interpolation of projection subspaces method. When employing these techniques, a trade-off should be found between accuracy and computation time.

### 3. EQS problem

#### 3.1. Finite element EQS problem

By neglecting the effect of the time variation of the magnetic flux density, the electro-quasistatic (EQS) problem in the frequency domain is defined by:

$$ \mathrm{curl} E = 0, \qquad (25) $$

$$ \mathrm{div}(\mathbf{J} + j\omega\mathbf{D}) = 0. \qquad (26) $$

with $E$ the electric field, $\mathbf{J}$ the current density, $\mathbf{D}$ the electric displacement field, and $\omega$ the angular frequency. In the whole domain $D$, the constitutive relationship between the electric displacement field and the electric field is given by $\mathbf{D} = \epsilon\mathbf{E}$. The current density and the electric field are also linked by $\mathbf{J} = \sigma\mathbf{E}$. In the EQS problem, the conductivity in some regions of $D$ is so small that it can be considered equal to zero. The electric scalar potential $\varphi$ is generally introduced to solve the problem and the electric field is then written under the form:

$$ \mathbf{E} = -\mathrm{grad}\varphi - V\mathrm{grad}\alpha \qquad (27) $$

where the source field $\alpha$ is introduced to impose the homogeneous boundary conditions on the scalar potential. Combining (26) and (27), the problem to solve becomes:

$$ \mathrm{div}(\sigma\mathrm{grad}\varphi + j\omega(\varepsilon\mathrm{grad}\varphi) + \sigma V\mathrm{grad}\alpha + j\omega(\varepsilon V\mathrm{grad}\alpha)) = 0 \qquad (28) $$

If the source is the current $I$ imposed by the external circuit, the voltage becomes then an unknown. An additional equation should be added to solve the problem.

$$ I = - \int_{D_{\varphi}} \mathrm{grad}\alpha \cdot \mathbf{J}d\mathbf{D} - \int_{D} \mathrm{grad}\alpha \cdot j\omega \mathbf{D}d\mathbf{D} \qquad (29) $$

When applying the Finite Element method, the scalar potential $\varphi$ is discretized with nodal elements. We denote by $\mathbf{X}_{\varphi}$ the vector of Degrees of Freedom (DoF) of $\varphi$ of size $N_{\varphi}$. These DoFs are the values of the scalar potential at the nodes of the mesh. We have then:

$$ \varphi = \sum_{i=1}^{N_{\varphi}} w_{i}^{0}(\mathbf{x})X_{\varphi i} \qquad (30) $$

with $X_{\varphi i}$ being the $i$th component of $\mathbf{X}_{\varphi}$ and $w_{i}^{0}(\mathbf{x})$ the nodal interpolation function associated to the node $i$. Applying the Galerkin method to (28) and (29) results in the following equation system:

$$ \mathbf{G}_{w^{0},u^{0},\alpha}^{\sigma} = \left[ \int_{D} \mathrm{grad} w_{i}^{0} \cdot \sigma \mathrm{grad} w_{j}^{0} d\mathbf{D} \right]_{1 \le i, j \le N_{\varphi}} = -V \mathbf{F}_{w^{0},u^{0},\alpha}^{\sigma} - j\omega V \mathbf{F}_{w^{0},u^{0},\alpha}^{\varepsilon}, \qquad (31) $$

$$ [\mathbf{F}_{w^{0},\alpha}^{\sigma}]^{\mathrm{T}} \mathbf{X}_{\varphi} + \mathbf{G}_{a,a}^{\sigma} V + j\omega([\mathbf{F}_{w^{0},\alpha}^{\varepsilon}]^{\mathrm{T}} \mathbf{X}_{\varphi} + \mathbf{S}_{a,a}^{\varepsilon} V) = I. \qquad (32) $$

With $I$ the input current and the entries of the different matrices given by:

$$ \left\{ \begin{array}{l} \mathbf{G}_{w^0, u^0, \alpha}^{\sigma} = \left[ \int_D \mathrm{grad} w_i^0 \cdot \sigma \mathrm{grad} w_j^0 d\mathbf{D} \right]_{1 \le i, j \le N_{\varphi}} \\ \mathbf{S}_{w^0, u^0, \alpha}^{\varepsilon} = \left[ \int_D \mathrm{grad} w_i^0 \cdot \varepsilon \mathrm{grad} w_j^0 d\mathbf{D} \right]_{1 \le i, j \le N_{\varphi}} \\ \mathbf{F}_{w^0, \alpha}^{\sigma} = \left[ \int_D \mathrm{grad} w_i^0 \cdot \sigma \mathrm{grad} \alpha d\Omega \right]_{1 \le i \le N_{\varphi}} \\ \mathbf{F}_{w^0, \alpha}^{\varepsilon} = \left[ \int_D \mathrm{grad} w_i^0 \cdot \varepsilon \mathrm{grad} \alpha d\Omega \right]_{1 \le i \le N_{\varphi}} \\ \mathbf{G}_{a,a}^{\sigma} = \int_D \mathrm{grad} \alpha \cdot \sigma \mathrm{grad} \alpha d\mathbf{D} \\ \mathbf{S}_{a,a}^{\varepsilon} = \int_D \mathrm{grad} \alpha \cdot \varepsilon \mathrm{grad} \alpha d\mathbf{D} \end{array} \right. \qquad (33) $$
<figure><img src="image_4.png" /><figcaption>Fig. 1. (a) Illustration (not to scale) of the symmetry plane in the simulation domain for a healthy bushing, highlighting the applied voltage source, the resulted current flow, and corresponding boundary conditions. (b) Location of the defect.</figcaption></figure>

Gathering the vector $\mathbf{X}_{\phi}$ and the voltage $V$ into one unknown vector:

$$ \mathbf{X}_{\phi V} = \begin{bmatrix} \mathbf{X}_{\phi} \\ V \end{bmatrix} \qquad (34) $$

The equation system (31) and (32) are gathered under the form:

$$ \mathbf{N}_{\phi V} \mathbf{X}_{\phi V} + j\omega \mathbf{K}_{\phi V} \mathbf{X}_{\phi V} = I \mathbf{F}_{\phi V}, \qquad (35) $$

with the matrices $\mathbf{N}_{\phi V}$ and $\mathbf{K}_{\phi V}$ and vector $\mathbf{F}_{\phi V}$ defined as follows:

$$ \mathbf{N}_{\phi V} = \begin{bmatrix} \mathbf{G}_{\omega}^{\sigma} & \mathbf{F}_{\omega}^{\sigma} \\ [\mathbf{F}_{\omega}^{\sigma}, \omega]^{T} & \mathbf{G}_{\omega, \omega}^{\sigma} \end{bmatrix}, \qquad (36) $$

$$ \mathbf{K}_{\phi V} = \begin{bmatrix} \mathbf{S}_{\omega}^{\epsilon} & \mathbf{F}_{\omega}^{\epsilon} \\ [\mathbf{F}_{\omega}^{\epsilon}, \omega]^{T} & \mathbf{S}_{\omega, \omega}^{\epsilon} \end{bmatrix}, \qquad (37) $$

$$ \mathbf{F}_{\phi V} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}. \qquad (38) $$

As mentioned previously, the conductivity $\sigma$ can be equal to zero in some regions leading to a singular matrix $\mathbf{N}_{\phi V}$. When applying model order reduction technic, [31] proposed to make a change of variable $s = \frac{1}{j\omega}$ leading to the final system:

$$ s\mathbf{N}_{\phi V}\mathbf{X}_{\phi V} + \mathbf{K}_{\phi V}\mathbf{X}_{\phi V} = sI\mathbf{F}_{\phi V}, \qquad (39) $$

We obtain the form similar to (1), with $\mathbf{K} = \mathbf{K}_{\phi V}$, $\mathbf{N} = \mathbf{N}_{\phi V}$, and $\mathbf{F} = \mathbf{F}_{\phi V}$.

## 3.2. Parametric EQS problem

We will now consider a parametric dependence considering the problem introduced in [31]. The studied domain holds a Resin-Impregnated Paper (RIP) bushing and the surrounding air. The RIP bushing is made of nine stratified insulation layers, each separated by a fine capacitor screen that we will consider equipotential. If there is no defect, the problem is axially symmetric as shown in Fig. 1(a). It can appear localized inter-screens insulation degradation and we want to study the impact of such a defect on the behavior of the bushing. A high voltage $V$ appears between the inner face and the outer face of the bushing, and a current $I$ appears that is rather small under healthy condition but increases when an insulation default arises. In order to be able to predict a failure, an indicator has been proposed to evaluate the level of insulator degradation which is defined as follows:

$$ \tan \delta(\omega) = -\frac{V_R}{V_C}, \quad \text{with } V = V_R + jV_C \qquad (40) $$

with $V$ the voltage when a current $I = 1$ A is imposed.

We assume that there is a defect in the eighth layer of the bushing. The defect is a cylinder with a diameter of 5.7 mm, which connects the capacitor screens of the eighth layer of the bushing as shown in Fig. 1(b). The symmetry of the problem is broken, which then becomes 3D.
The boundary conditions are shown in Fig. 1(a). A current I is injected through the inner face and it appears then a voltage V between the inner face and the outer face of the bushing.

The material behavior law is assumed to be frequency-invariant and linear. The electrical conductivity and relative permittivity of the bushing are denoted respectively $\sigma_n$ and $\epsilon_n$, respectively. The conductivity is assumed to be constant and equal to $\sigma_n = 8.33 \times 10^{-13}$ S/m. The permittivity $\epsilon_n$ varies with the humidity in the RIP. Hence, it is considered as a parameter for the parametric problem and its relative value varies from 3.156 to 20.157. The conductivity of the air surrounding the bushing is assumed to be 0 and the relative permittivity is equal to 1. The relative permittivity of the defect is equal to 11.66. The conductivity of the defect is assumed to be equal to $10^5 \cdot \sigma_n$ S/m. In practice, we do not know the exact value of the conductivity in the defect and we also want to study the influence of this parameter. We assume that the variation can be very high and that the parameter $\zeta$ can vary from $\zeta_{min} = 1$ to $\zeta_{max} = 9$. We have then finally two parameters the permittivity of the RIP and the conductivity of the defect.

The parameterized full FE system of this application is then obtained by incorporating the dependence of the parameter $\zeta$ and $\epsilon_n$ into the equation system (39). We have then:

$$s\mathbf{N}_{\phi V}(\zeta)\mathbf{X}_{\phi V}(s, \zeta, \epsilon_n) + \mathbf{K}_{\phi V}(\epsilon_n)\mathbf{X}_{\phi V}(s, \zeta, \epsilon_n) = s\mathbf{I}\mathbf{F}_{\phi V} \quad (41)$$

The matrix $\mathbf{N}_{\phi V}(\zeta)$ is a function of the conductivity, so a function of the parameter $\zeta$. Similarly, the matrix $\mathbf{K}_{\phi V}(\epsilon_n)$ is a function of the parameter $\epsilon_n$. We want to apply now the two previous methods, the MultiPMomMtch method and the projection subspaces method.

### 3.3. Multi-moment matching method

#### 3.3.1. Affine decomposition of the EQS problem

To apply the multi-moment matching method, the parametric problem should be written under an affine decomposition. Since the parameters are all related to the material parameters, the affine decomposition of the matrix $\mathbf{N}_{\phi V}(\zeta)$ and $\mathbf{K}_{\phi V}(\epsilon_n)$ can be obtained by splitting the integration domain into different parts. In particular, we have the affine decomposition for $\mathbf{N}_{\phi V}(\zeta)$ which can be written under the form:

$$\mathbf{N}_{\phi V}(\zeta) = \mathbf{N}_{\phi V,0} + 10^{\zeta} \cdot \mathbf{N}_{\phi V,1} \quad (42)$$

The expressions of the matrices $\mathbf{N}_{\phi V}^{\sigma_n}$ and $\mathbf{N}_{\phi V}^{\sigma_d}$ are:

$$\mathbf{N}_{\phi V,0} = \begin{bmatrix} \left[ \int_{D_n} \mathbf{grad} w_i(x)^0 \cdot \sigma(x) \mathbf{grad} w_j^0(x) d\Omega \right]_{1 \le i \le N_{\phi}} & \left[ \int_{D_n} \mathbf{grad} w_i^0 \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \right]_{1 \le i \le N_{\phi}} \\ \left[ \int_{D_n} \mathbf{grad} w_i^0 \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \right]^{T}_{1 \le i \le N_{\phi}} & \int_{D_n} \mathbf{grad} \alpha \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \end{bmatrix}, \quad (43)$$

$$\mathbf{N}_{\phi V,1} = \begin{bmatrix} \left[ \int_{D_d} \mathbf{grad} w_i(x)^0 \cdot \sigma(x) \mathbf{grad} w_j^0(x) d\Omega \right]_{1 \le i \le N_{\phi}} & \left[ \int_{D_d} \mathbf{grad} w_i^0 \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \right]_{1 \le i \le N_{\phi}} \\ \left[ \int_{D_d} \mathbf{grad} w_i^0 \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \right]^{T}_{1 \le i \le N_{\phi}} & \int_{D_d} \mathbf{grad} \alpha \cdot \sigma(x) \mathbf{grad} \alpha d\Omega \end{bmatrix}. \quad (44)$$

Where $D_d$ denotes the domain of the defect and $D_n$ is defined as $D_n = D \setminus D_d$ knowing that D is the whole computational domain. The conductivity $\sigma(x)$ in (43) and (44) is equal to $10^5 \cdot \sigma_n$ in the defect, $\sigma_n$ in the RIP, and zero elsewhere. We can see that the matrices $\mathbf{N}_{\phi V,0}$ and $\mathbf{N}_{\phi V,1}$ are similar to $\mathbf{N}_{\phi V}$ except that they are constructed only in the domains $D_n$ and $D_d$, respectively.

We have also an affine decomposition for $\mathbf{K}_{\phi V}(\epsilon_n)$ which can be written under the form:

$$\mathbf{K}_{\phi V}(\epsilon_n) = \mathbf{K}_{\phi V,0} + \mathbf{K}_{\phi V,1}(\epsilon_n) \quad (45)$$

The expression of matrices $\mathbf{K}_{\phi V,0}$ and $\mathbf{K}_{\phi V,1}(\epsilon_n)$ are:

$$\mathbf{K}_{\phi V,0} = \begin{bmatrix} \mathbf{S}_{w^{0,0}}^{e} & \mathbf{F}_{0}^{e} \\ [\mathbf{F}_{0}^{e}]^{T} & \mathbf{S}_{\alpha,0}^{e} \end{bmatrix}, \quad (46)$$

$$\mathbf{K}_{\phi V,1}(\epsilon_n) = \begin{bmatrix} \mathbf{S}_{w^{0,1}}^{e}(\epsilon_n) & \mathbf{F}_{1}^{e}(\epsilon_n) \\ [\mathbf{F}_{1}^{e}(\epsilon_n)]^{T} & \mathbf{S}_{\alpha,1}^{e}(\epsilon_n) \end{bmatrix}. \quad (47)$$

With the entries of the different matrices given by:

$$\left\{ \begin{array}{l} \mathbf{S}_{w^{0,0}}^{e} = \left[ \int_{D_s} \mathbf{grad} w_i(x)^0 \cdot \epsilon(x) \mathbf{grad} w_j^0(x) d\Omega \right]_{1 \le i \le N_{\phi}} \\ \mathbf{F}_{0}^{e} = \left[ \int_{D_s} \mathbf{grad} w_i^0 \cdot \epsilon(x) \mathbf{grad} \alpha d\Omega \right]_{1 \le i \le N_{\phi}} \\ \mathbf{S}_{\alpha,0}^{e} = \int_{D_s} \mathbf{grad} \alpha \cdot \epsilon(x) \mathbf{grad} \alpha d\Omega \\ \mathbf{S}_{w^{0,1}}^{e}(\epsilon_n) = \left[ \int_{D_{RIP}} \mathbf{grad} w_i(x)^0 \cdot \epsilon(x) \mathbf{grad} w_j^0(x) d\Omega \right]_{1 \le i \le N_{\phi}} \\ \mathbf{F}_{1}^{e}(\epsilon_n) = \left[ \int_{D_{RIP}} \mathbf{grad} w_i^0 \cdot \epsilon(x) \mathbf{grad} \alpha d\Omega \right]_{1 \le i \le N_{\phi}} \\ \mathbf{S}_{\alpha,1}^{e}(\epsilon_n) = \int_{D_{RIP}} \mathbf{grad} \alpha \cdot \epsilon(x) \mathbf{grad} \alpha d\Omega \end{array} \right. \quad (48)$$
<table><caption>**Table 2**
Dimension of reduced bases obtained by the MultiPMomMtch method and the interpolation of projection subspaces method.</caption><thead><tr><th></th><th>Size</th><th>Number of bases</th></tr></thead><tbody><tr><td>Base V for the MultiPMomMtch method</td><td>9266 × 30</td><td>1</td></tr><tr><td>Base W(p) Interpolation of projection subspaces</td><td>9266 × 10</td><td>6</td></tr></tbody></table>

where $D_{RIP}$ denotes the domain of the RIP excluding the defect and $D_s = D \setminus D_{RIP}$. The permittivity $\epsilon(\mathbf{x})$ in (46) and (47) is equal to $\epsilon_n$ in the RIP, 11.66$\epsilon_0$ in the defect, and $\epsilon_0$ elsewhere where $\epsilon_0$ is the permittivity of the vacuum.

We denote $\hat{\mathbf{p}} = (s, \zeta, \epsilon_n)$. From (42) and (45), we can deduce an affine decomposition which is required if we want to apply the MultiPMomMtch method:

$$ \mathbf{Q}(\hat{\mathbf{p}}) = \mathbf{Q}(s, \zeta, \epsilon_n) = \mathbf{K}_{\phi V,0} + \mathbf{K}_{\phi V,1}(\epsilon_n) + s\mathbf{N}_{\phi V,0} + s \cdot 10^{\zeta} \cdot \mathbf{N}_{\phi V,1} \quad (49) $$

### 3.3.2. Configuration for the MultiPMomMtch method in the offline stage

The $l = 8$ expansion points are chosen as the vertices of the hypercube formed by the domain of variation of the two parameters and the parameter $s$:

$$ \begin{gathered} \zeta \in [1, 9], \\ \epsilon_n \in [3.156, 20.157], \\ s \in [\frac{1}{2\pi} \times 10^{-3}, \frac{1}{2\pi} \times 10^{3}]. \end{gathered} \quad (50) $$

The order of the series expansion $q$ is equal to 1. In summary, the characteristic values (see Table 1 for their definitions) of the MultiPMomMtch method for the EQS problem are:

$$ d = 1, \quad p = 3, \quad q = 1, \quad l = 8. \quad (51) $$

We have applied $l = 2^3 = 8$ times the algorithm 1. The reduced basis V then has $(p+1)l = 32$ vectors. Since the reduced basis does not depend on the parameters $\mathbf{p}$, in the online stage, there is no additional cost to calculate the reduced basis.

## 3.4. Interpolation of projection subspaces method

### 3.4.1. Offline stage

In the offline stage, we have chosen to precompute $K = 6$ local basis matrices $\mathbf{W}_1, \mathbf{W}_2, \dots, \mathbf{W}_6$, each obtained by applying 5 times the full EQS problem. The sample of parameter $\mathbf{p} = (\zeta, \epsilon_n)$ has been chosen as 4 corners of the cube: $[\zeta_{min}, \zeta_{max}] \times [\epsilon_{min}, \epsilon_{max}]$ and another two midpoints in the cube $(\frac{\zeta_{max}+\zeta_{min}}{2}, \epsilon_{min})$ and $(\frac{\zeta_{max}+\zeta_{min}}{2}, \epsilon_{max})$. For example, the matrix $\mathbf{W}_1$ is constructed as:

$$ \mathbf{W}_1 = [\mathrm{Re}(\mathbf{X}_1) \ \mathrm{Re}(\mathbf{X}_2) \ \dots \ \mathrm{Re}(\mathbf{X}_5) \ \mathrm{Im}(\mathbf{X}_1) \ \mathrm{Im}(\mathbf{X}_2) \ \dots \ \mathrm{Im}(\mathbf{X}_5)] \quad (52) $$

where $\mathbf{X}_1, \dots, \mathbf{X}_5$ are solutions to the EQS problem (41) at 5 different frequency points equally spaced in the interval $[1 \text{ mHz}, 1 \text{ Khz}]$ in log scale.

### 3.4.2. Online stage

In the online stage, the weight function used is the multi linear interpolation function. The contact point is the subspace spanned by the basis which is constructed at parameter $\mathbf{p}_k$ and $\mathbf{p}_k$ is the nearest parameter point to the parameter $\mathbf{p}$.

Table 2 shows the size of reduced bases obtained by the two PMOR methods. We can see that the number of full problems solved for the two methods during the offline stage is the same by counting the number of columns of the reduced bases. That is 30 for the MultiPMomMtch method, and 60/2 for the interpolation of projection subspaces method since each basis is composed of real part and imaginary part of the solutions.

## 3.5. Comparison

Fig. 2 shows the variation of the indicator $\tan(\delta)$ calculated with the full EQS problem in function of the frequency $f$ and the parameter $\zeta$. We can see a strong variation of this indicator with the frequency as well as a non-monotonous variation with the parameter $\zeta$.

In the following, we will compare the full EQS model with the two reduced models obtained by the MultiPMomMtch method and the interpolation of projection subspaces method at $R = 36$ points of the parameter set $\mathbf{p}$. These $R$ points are randomly distributed and generated by applying the Latin Hypercube method. The full model will be taken as the reference. The errors of the indicator between the full model and the two reduction methods are given by:

$$ \Delta_{\text{rel,tan}(\delta)}(\mathbf{p}) = \max_{1 \le i \le 20} \frac{|\tan(\delta)(\omega_i) - \tan(\delta)(\omega_i)|}{\tan(\delta)(\omega_i)} \quad (53) $$
<figure><img src="image_5.png" /><figcaption>Fig. 2. Variation of tan(δ) obtained by the EQS problem with ζ, where f is the frequency.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Histogram of Δrel,tan(δ)(p) obtained the MultiPMomMtch method and the interpolation of projection subspaces method, where 36 Latin hypercube samples of p are chosen.</figcaption></figure>

where the angular frequencies ωᵢ are logarithmically distributed in the interval [2π×10⁻³, 2π×10³], tan(δ)(ωᵢ) are the results obtained by the PROM methods at ωᵢ and p.

Similarly, we can calculate the error of the local solution between the full model and the two reduction methods, as shown in Fig. 4:

$$ \Delta_{\text{rel},\mathbf{x}}(\mathbf{p}) = \max_{1 \le i \le 20} \frac{\|\mathbf{X}_i - \hat{\mathbf{X}}_i\|_2}{\|\mathbf{X}_i\|_2} (\mathbf{p}) \quad (54) $$

where $\hat{\mathbf{X}}_i$ is the approximated local solution obtained by the PROM method at ωᵢ and p and $\mathbf{X}_i$ is the local solution at ωᵢ and p.

In Fig. 3, we give the histogram of the log of the error Δrel,tan(δ)(p) on the global quantity tan(δ) obtained from R = 36 value of errors calculated from (53). It appears clearly that the MultiMomMtch method is the most accurate. In Fig. 4, we present a comparison of the log of the error obtained on the solution of the problem given by (54). Again, the MultiMomMtch method remains the most accurate. Compared to the interpolation of projection subspaces method, the multi-moment matching method is more accurate because it takes advantage of the affine decomposition, which enables a Taylor series approximation of the transfer function around the expansion point

Figs. 5 and 6 show the variation of Δrel,tan(δ)(p) with ζ and εₙ, respectively, when the other parameter is fixed. Since the interpolation of projection subspaces method is interpolant, and if there is no truncation when applying the SVD on the snapshot in the offline stage, then if the parameter set is equal to one of the parameter sets used during the offline stage for the snapshot determination, the error is equal to zero. If a parameter p generated randomly is close to a parameter used to during the offline stage, then the error will be small. This is why we can see that, when ζ is close to 1 or 9 and when εₙ is close to 3.156 or 20.157,
<figure><img src="image_7.png" /><figcaption>Fig. 4. Histogram of $\Delta_{rel,x}(\mathbf{p})$ obtained the MultiPMomMtch method and the interpolation of projection subspaces method.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Variation of $\Delta_{rel,tan(\delta)}(\mathbf{p})$ obtained by the MultiPMomMtch method (on the left) and the interpolation of projection subspaces method (on the right), with $\epsilon_n$ fixed at 13.02.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Variation of $\Delta_{rel,tan(\delta)}(\mathbf{p})$ obtained by the MultiPMomMtch method (on the left) and the interpolation of projection subspaces method (on the right), with $\zeta$ fixed at 7.5.</figcaption></figure>

this method is the most accurate. It shows also that this method is influenced by the *K* parameter points chosen to calculate the snapshots in the offline stage.
<table><caption>**Table 3**

Speedup of the MultiPMomMtch method, and the interpolation of projection subspaces method for the EQS problem.</caption><thead><tr><th></th><th>Offline stage</th><th>Online stage</th></tr></thead><tbody><tr><td>MultiPMomMtch</td><td>1.3</td><td>2.1 × 10<sup>4</sup></td></tr><tr><td>Interpolation of projection subspaces</td><td>4.0</td><td>52</td></tr></tbody></table>

Table 3 presents the speedup for:

1. The construction of global basis matrices for the MultiPMomMtch method in the offline stage.

2. The reduced system to calculate 20 points in the frequency domain for the MultiPMomMtch method in the online stage under case $\mathbf{p} = (\zeta = 7.0453, \epsilon_n = 14.4375)$.

3. The construction of $K = 6$ local basis matrices for the interpolation of projection subspaces method in the offline stage.

4. The construction of the matrix $V$ for the interpolation of projection subspaces method in the online stage under case $\mathbf{p} = (\zeta = 7.0453, \epsilon_n = 14.4375)$.

Where the speedup is calculated by:

$$ \text{speedup} = \frac{T_{EF}}{T} \qquad (55) $$

with $T_{EF}$ the computational time for the full FE model to calculate 20 points in the frequency domain regularly distributed in the interval $[1 \times 10^{-3}, 1 \times 10^3]$ Hz in log scale.

In the online stage, the interpolation of projection subspaces method takes more time compared to the MultiPMomMtch method since 4 basis matrices need to be loaded and 4 SVD operations need to be carried out each time the parameter changes. On the other hand, for the MultiPMomMtch method, only the reduced system needs to be solved since the reduced basis is not parameter dependent.

# 4. MQS problem

## 4.1. Finite element MQS problem

The magneto-quasistatic (MQS) problem in the frequency domain is defined by the following equations:

$$ \text{curlH} = \mathbf{J}, \qquad (56) $$

$$ \text{curlE} = -j\omega\mathbf{B}, \qquad (57) $$

$$ \text{div} \mathbf{B} = 0, \qquad (58) $$

$$ \text{div} \mathbf{J} = 0. \qquad (59) $$

with **H** the magnetic field and **B** the magnetic flux density. The constitutive relationship between **H** and **B** is assumed to be linear, and we have **B** = $\mu\mathbf{H}$ where $\mu$ is the permeability. The second constitutive relationship is defined by the Ohm law $\mathbf{J} = \sigma\mathbf{E}$ with $\sigma$ the conductivity. The permeability is a function always strictly positive, whereas the conductivity can be equal to zero in nonconductive region. As we shall see later, the problem treated in this paper is composed of a thin structure. Solving this kind of problem for the potential formulation $\mathbf{A} - \varphi$ with the iterative solver is difficult if the structure is discretized with the standard three-dimensional element [32]. However, for the nonconductive domain, the matrix obtained from the vector potential $\mathbf{A}$ can be easily parameterized with the geometry parameter in our problem. To solve the problem, we will use the double vector potential formulation $\mathbf{A} - \mathbf{T}$ [33], which combines both the advantages of the formulation $\mathbf{A} - \varphi$ and the formulation $\mathbf{T} - \Omega$. We can express then the current density $\mathbf{J}$ and the magnetic flux density $\mathbf{B}$ as:

$$ \begin{aligned} \mathbf{J} &= I\lambda_I + \text{curlT}, \quad \text{in } D_c, \\ \mathbf{B} &= \text{curlA}, \quad \text{in } D. \end{aligned} \qquad (60) $$

Where $D_c$ denotes the conductive domain where the conductivity $\sigma$ is strictly positive, $D$ denotes the whole domain, and $\lambda_I$ is the source field [34] introduced to impose the source of current $I$ flowing through the conductor. We denote $V$ the voltage imposed by the external circuit to obtain the current $I$. The current is then a function of the frequency $\omega$.

Combining with (56) and (57), the strong problem is obtained:

$$ \text{curl} \frac{1}{\mu} \text{curlA} = \text{curlT} + I(\omega)\lambda_I, \qquad (61) $$

$$ \text{curl} \frac{1}{\sigma} (\text{curlT} + I(\omega)\lambda_I) = -j\omega \text{curlA}. \qquad (62) $$

The vector potentials $\mathbf{A}$ and $\mathbf{T}$ are discretized with edge elements in the whole domain and in the conductive domain, respectively. The source term $\lambda_I$ is expressed with facet elements. We denote by $\mathbf{X}_A$ the vector of size $N_A$ associated with the Dof of $\mathbf{A}$. Similarly,
<figure><img src="image_10.png" /><figcaption>Fig. 7. (a) Assembling of the planar coil and the ferrite core. (b) Planar coil only. (c) Depiction of the coil's winding orientation.</figcaption></figure>

we denote by $X_T$ the vector of size $N_T$ associated with the Dof's of $T$. The number of Dof's of $A$ and $T$ is about 2.0 million and 40 thousand, respectively. We have then:

$$ 
\begin{aligned}
\mathbf{A} &= \sum_{i=1}^{N_A} \mathbf{w}_i^1(\mathbf{x}) X_{A_i}, \\
\mathbf{T} &= \sum_{i=1}^{N_T} \mathbf{w}_i^1(\mathbf{x}) X_{T_i}.
\end{aligned} 
\qquad (63) $$

with $X_{A_i}$ is $i$th component of $X_A$ which corresponds to the circulation of $A$ along the edge $i$ and $w_i^1(\mathbf{x})$ its edge interpolation function. Similarly, we denote $X_{T_i}$ the $i$th component of $X_T$ which corresponds to the circulation of $T$ along the edge $i$. Applying the FE method to (56) and (57) results in (taking $s = j\omega$):

$$ \begin{bmatrix} \mathbf{K}_A & -\mathbf{L} \\ 0 & \mathbf{K}_T \end{bmatrix} \begin{bmatrix} \mathbf{X}_A \\ \mathbf{X}_T \end{bmatrix} + s \begin{bmatrix} 0 & 0 \\ \mathbf{L}^T & 0 \end{bmatrix} \begin{bmatrix} \mathbf{X}_A \\ \mathbf{X}_T \end{bmatrix} = I(s) \begin{bmatrix} \mathbf{F}_A \\ \mathbf{F}_T \end{bmatrix} \qquad (64) $$

where the entries of the matrices are given by:

$$ 
\left\{ 
\begin{array}{l}
\mathbf{L} = \left[ \int_{D_c} \mathbf{w}_i^1 \mathrm{curl} \mathbf{w}_j^1 dD \right]_{1 \le i \le N_A; 1 \le j \le N_T} \\
\mathbf{F}_A = \left[ \int_{D_c} \mathbf{w}_i^1 \cdot \lambda_I dD \right]_{1 \le i \le N_A} \\
\mathbf{F}_T = \left[ -\int_D \mathbf{w}_i^1 \cdot \frac{1}{\sigma} \lambda_I dD \right]_{1 \le i \le N_T}
\end{array}
\right. 
\qquad (65) $$

Equation system (64) is not symmetric. In the following, we will show how to obtain a symmetric one. First, an electrokinetic problem can be solved:

$$ \mathbf{K}_T \mathbf{X}_T = \mathbf{F}_T \qquad (66) $$

Let the solution of (66) be $\mathbf{X}_{T,0}$ which can be used for the change of variable:

$$ \mathbf{X}_T = \mathbf{X}_{T,1} + I(s)\mathbf{X}_{T,0} \qquad (67) $$

where $\mathbf{X}_{T,1}$ is unknown. Then Eq. (64) can be reformulated as:

$$ \begin{bmatrix} \mathbf{K}_A & -\mathbf{L} \\ 0 & \mathbf{K}_T \end{bmatrix} \begin{bmatrix} \mathbf{X}_A \\ \mathbf{X}_{T,1} \end{bmatrix} + s \begin{bmatrix} 0 & 0 \\ \mathbf{L}^T & 0 \end{bmatrix} \begin{bmatrix} \mathbf{X}_A \\ \mathbf{X}_{T,1} \end{bmatrix} = I(s) \begin{bmatrix} \mathbf{b}_{AT} \\ \mathbf{0} \end{bmatrix} \qquad (68) $$

where the vector $\mathbf{b}_{AT}$ is defined as:

$$ \mathbf{b}_{AT} = \mathbf{F}_A + \mathbf{L}\mathbf{X}_{T,0} \qquad (69) $$

By replacing the vector $\mathbf{X}_{T,1}$ by its expression in function of the vector $\mathbf{X}_A$, we can obtain from (68):

$$ \mathbf{K}_A \mathbf{X}_A + s \mathbf{L} \mathbf{K}_T^{-1} \mathbf{L}^T \mathbf{X}_A = \mathbf{b}_{AT} I(s) \qquad (70) $$

The quantity of interest for this problem is the voltage $V$ and we have:

$$ V(s) = s \mathbf{b}_{AT}^T \mathbf{X}_A + R_0 I(s) \qquad (71) $$

where the constant $R_0$ and the quantity $R_0 I$ is the effect associated with the DC current. Once again, a form similar to (1) is obtained, with $\mathbf{K} = \mathbf{K}_A$, $\mathbf{N} = \mathbf{L} \mathbf{K}_T^{-1} \mathbf{L}^T$, and $\mathbf{F} = \mathbf{b}_{AT}$. The global quantity that will considered in the following is the impedance give by:

$$ Z(s) = \frac{V(s)}{I(s)} = R(s) + sL(s) \qquad (72) $$

where the function $R(s)$ represents the resistance and the function $L(s)$ is the inductance.
<figure><img src="image_11.png" /><figcaption>Fig. 8. Partitions of the simulated domain.</figcaption></figure>

## 4.2. Parametric MQS problem

We consider the planar inductor presented in Fig. 7. This inductor is a component of power electronic converters. During the design stage of such planar inductor, the influence of parameters like permeability of the iron core, the size of the airgap ..., should be studied in order to optimize the converter, then this optimization requires the simulation of the planar inductor with the external circuit. This simulation cannot be done directly using a Finite Element model because it would be too time consuming. This is why a reduced model should be developed. The parameters of the planar inductor, we consider in the following, are the relative permeability of the upper ferrite core $\mu_1$, the thickness of the air gap $\delta_1$ and the position of the conductor in the window $\delta_2$ (see Fig. 8). The air gap is meshed with 5 layers of prism elements extruded from a surface meshed with triangular elements. If the parameter $\delta_1$ changes, the mesh of the air gap will also change by modifying the shape of the elements. If we incorporate the parameter dependence into the full system (70), we can obtain the parameterized system:

$$ \mathbf{K}_A(\delta_1, \delta_2, \mu_1)\mathbf{X}_A(\delta_1, \delta_2, \mu_1, s) + s\mathbf{L}\mathbf{K}_T^{-1}\mathbf{L}^T\mathbf{X}_A(\delta_1, \delta_2, \mu_1, s) = \mathbf{b}_{AT}I(s) \quad (73) $$

### 4.2.1. MultiPMomMtch method

To apply the MultiPMomMtch method, we need an affine decomposition of the matrix $\mathbf{K}_A(\delta_1, \delta_2, \mu_1)$. We can show in Appendix B that only the matrix $\mathbf{K}_A(\delta_1, \delta_2, \mu_1)$ is a function of $\delta_1$, $\delta_2$ and $\mu_1$, and that we still have an affine decomposition:

$$ \mathbf{K}_A(\delta_1, \delta_2, \mu_1) = \mathbf{K}_0 + \frac{1}{\delta_1} \cdot \mathbf{K}_1 + \delta_1 \cdot \mathbf{K}_2 + \frac{1}{\delta_2} \cdot \mathbf{K}_3 + \delta_2 \cdot \mathbf{K}_4 + \frac{1}{2\delta_1 - \delta_2} \cdot \mathbf{K}_5 + (2\delta_1 + \delta_2) \cdot \mathbf{K}_6 + \mu_1 \mathbf{K}_7 \quad (74) $$

The order of series expansion $q$ is equal to 1. In summary, the characteristic values (see Table 1 for their definition) of the MultiPMomMtch method for the MQS problem are:

$$ d = 1, \quad p = 8, \quad q = 1, \quad l = 16. \quad (75) $$

The expansion points are chosen as the vertices of the hypercube formed by the domain of variation of the three parameters defined previously and the parameter $s$:

$$ \begin{array}{l} \delta_1 \in [0.1, 1], \\ \delta_2 \in [0.1, 2.5], \\ \mu_1 \in [200, 2000], \\ s \in [2\pi \times 10^3, 10\pi \times 10^5]. \end{array} \quad (76) $$

We have applied $l = 2^4 = 16$ times Algorithm 1. The global basis matrix $V$ then has $(p+1)l = 144$ vectors.

## 4.3. Interpolation of projection subspaces method

In the offline stage, we need to calculate the snapshots. To span the parameter space, we have considered the following values for the three parameters:

$$ \begin{array}{l} \delta_1 \in \{0.1, 0.55, 1\}, \\ \delta_2 \in \{0.1, 1.3, 2.5\}, \\ \mu_1 \in \{200, 2000\}. \end{array} \quad (77) $$

It leads to $K = 18$ possible sets of parameters, for which we have computed $K = 18$ local basis matrices $\mathbf{W}_1, \mathbf{W}_2, \dots, \mathbf{W}_{18}$, each basis $\mathbf{W}_i$ is constructed by solving the full MQS problem for 8 different frequencies to obtain the solution vectors: $\mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_8$. Specifically, $\mathbf{W}_i$ is formed by columns which are neither the real part nor the imaginary part of the solution vectors:

$$ \mathbf{W}_i = [\text{Re}(\mathbf{X}_1) \quad \dots \quad \text{Re}(\mathbf{X}_8) \quad \text{Im}(\mathbf{X}_1) \quad \dots \quad \text{Im}(\mathbf{X}_8)] \quad (78) $$
<table><caption>**Table 4**
Dimension of reduced bases obtained by the MultiPMomMtch method and the interpolation of projection subspaces method.</caption><thead><tr><td></td><td>Size</td><td>Number of bases</td></tr></thead><tbody><tr><td>Base V for the MultiPMomMtch method</td><td>1 952 932 × 144</td><td>1</td></tr><tr><td>Base W(p) Interpolation of projection subspaces</td><td>1 952 932 × 16</td><td>18</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Log of error of the resistance (left) and the inductance (right) for the MultiPMomMtch method, the interpolation of projection subspaces method without structure preserving and the interpolation of projection subspaces method with structure preserving.</figcaption></figure>

In total, we need to solve 18 × 8 = 144 times the full MQS problems which is the same number of MQS problems needed to be solved for the MultiPMomMtch method.

Table 4 shows the size of reduced bases obtained by the two PMOR methods. Like for the EQS application, we can see that the number of full problems solved for the two methods is the same.

### 4.4. Comparison

To compare the two methods of calculation of reduced basis, we will consider the frequency response of the impedance of the planar coil, and more specifically, the resistance and the inductance (see (72)). We will consider the FE MQS problem as a reference and calculate the error between the resistances and inductances obtained with the MultiPMomMtch method and the interpolation of projection subspaces method at 36 values of the parameter set ($\delta_1, \delta_2, \mu$) chosen randomly. Those values are generated using the Latin hypercube sampling method. For each point of the parameter set, we will calculate the resistance and inductance at 6 different frequencies $f \in [1 \times 10^3, 5 \times 10^5]$ to calculate the maximum relative error of resistance and inductance between the full MQS model and the PMOR model.

The comparison of the interpolation of projection subspaces method and the MultiPMomMtch method is shown in Fig. 9. We can see that the interpolation of projection subspaces method is less accurate than the MultiPMomMtch method for the 36 points with the definition of the error on the resistance and the inductance:

$$ \Delta_{\text{rel, R}}(\mathbf{p}) = \max_{1 \le i \le 6} \frac{|R(\omega_i) - \hat{R}(\omega_i)|}{R(\omega_i)}(\mathbf{p}), \qquad (79) $$

$$ \Delta_{\text{rel, L}}(\mathbf{p}) = \max_{1 \le i \le 6} \frac{|L(\omega_i) - \hat{L}(\omega_i)|}{L(\omega_i)}(\mathbf{p}). \qquad (80) $$

Where the angular frequencies $\omega_i$ are logarithmically distributed in the interval $[2\pi \times 10^3, \pi \times 10^6]$, $\hat{R}(\omega_i)$ and $\hat{L}(\omega_i)$ are the values of the resistances and inductances obtained by the PROM method at $\omega_i$ and $\mathbf{p}$.
<table><caption>**Table 5**

Speedup of the MultiPMomMtch method, the interpolation of projection subspaces method without structure preserving and the interpolation of projection subspaces method with structure preserving for the MQS problem.</caption><thead><tr><th></th><th>Offline stage</th><th>Online stage</th></tr></thead><tbody><tr><td>MultiPMomMtch</td><td>5.4 × 10<sup>-2</sup></td><td>3.6 × 10<sup>6</sup></td></tr><tr><td>Interpolation of projection subspaces without structure preserving</td><td>5.9 × 10<sup>-2</sup></td><td>6.2 × 10<sup>2</sup></td></tr><tr><td>Interpolation of projection subspaces with structure preserving</td><td>5.9 × 10<sup>-2</sup></td><td>6.3 × 10<sup>2</sup></td></tr></tbody></table>

It is well known that when applying a model order reduction model, the structure of the full problem should be preserved. It means that in the MQS case, we have two kinds of DofS: those related to the vector potential T and those related to the vector potential A. We imagine increasing the precision of the interpolation of projection subspaces method by first splitting the reduced basis $W_i$ with $i = \{1, 2, ..., K\}$ into two parts:

$$W_i = \begin{bmatrix} Wt_i \\ Wa_i \end{bmatrix} \qquad (81)$$

where the matrix $Wt_i$ is of size $N_T \times 16$ and the matrix $Wa_i$ is of size $N_A \times 16$. The integer $N_T$ is the number of unknowns for the vector potential T and the integer $N_A$ is the number of unknowns for the vector potential A. Then, the matrices $Wt_1, ..., Wt_K$ are interpolated using the interpolation of projection subspaces method to get a new reduced basis $Wt(p)$. Similarly, the matrices $Wa_1, ..., Wa_K$ are interpolated using the interpolation of projection subspaces method to get a new reduced basis $Wa(p)$. The final projection matrix is obtained by regrouping the two matrices:

$$W(p) = \begin{bmatrix} Wt_i(p) \\ Wa_i(p) \end{bmatrix} \qquad (82)$$

The comparison of the full model and the interpolation of projection subspaces method with splitting bases is shown in Fig. 9. We can see that the structure preserving technique enables to improve the accuracy for the resistance and the inductance. However, the MultiPMomMtch method remains always the most accurate.

Table 5 presents the speedup for:

1. The construction of global basis matrices for the MultiPMomMtch method in the offline stage.

2. The reduced system to calculate 6 points in the frequency domain for the MultiPMomMtch method in the online stage under case $p = (0.4638, 1.135, 1929)$.

3. The construction of $K = 18$ local basis matrices for the interpolation of projection subspaces method in the offline stage.

4. The construction of the matrix $V$ for the interpolation of projection subspaces method in the online stage under case $p = (0.4638, 1.135, 1929)$.

Where $T_{EF}$ the computational time for the full FE model to calculate 6 points in the frequency domain regularly distributed in the interval [1, 500] kHz in log scale.

In the online stage, for the interpolation of projection subspaces method, 8 basis matrices need to be loaded and 8 SVD operations need to be carried out each time the parameter changes. We can see that compared to the EQS problem, the speedup of the interpolation of projection subspaces method decreases due to the larger number of Dof. Using the structure-preserving technique, the computational cost almost does not change since the Dof of vector potential T is much smaller than the one of the vector potential A. In the offline stage, the speedup of the MultiPMomMtch method is quite similar to the speedup of the interpolation of projection subspaces method since the same number of full FE problems need to be solved. In Table 5, the time of the offline stage of the MultiPMomMtch method is 20 times longer than the time $T_{FE}$ to calculate the frequency response with the FE model. We can see also that the online stage is 3.6 $10^6$ faster so the time computation can be considered negligible. So, if the reduced model constructed by the MultiPMomMtch method is used more that 20 times during an online stage, like a design process for example, it is worth using the MultiPMomMtch method. This situation is often the case in practice, especially when the design concerns the whole power converter meaning that other components, like power switches and other passive components, are also considered.

The distribution of the current density modulus at 500 kHz obtained by the MQS problem is shown in Fig. 10 for:

1. $\delta_2 = 2.5$, the PCB is located at the bottom of the window of the iron core,

2. $\delta_2 = 2.5$, the PCB is located at the top of the window close to the air gap.

We can see that as the PCB moves closer to the air gap (decrease of $\delta_2$), the current density modulus increases on the inner PCB border due to the magnetic fringing effect meet close to the airgap.

The distribution of the current density has been also calculated with the MultiPMomMtch and the interpolation of projection subspaces methods. The error distribution of the current density obtained with the MultiPMomMtch method and the MQS problem is shown in Fig. 11. Similarly, to the one between the interpolation of projection subspaces method with structure preserving and the MQS problem is shown in Fig. 12. We can see that the MultiPMomMtch method approximates better the local quantities than the interpolation of projection subspaces method which is in coherence with the previous results obtained for the global quantities, i.e. the resistance and the inductance.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Distribution of the current density modulus (A/m²) obtained by the MQS problem (|J<sub>MQS</sub>|) at 500 kHz and source current of 1 A for δ₁ = 0.2, δ₂ = 2.5, and μ₁ = 2000 (left) and for δ₁ = 0.2, δ₂ = 0.1, and μ₁ = 2000 (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Distribution of error (|J<sub>PMOR</sub> - J<sub>MQS</sub>|) on current density obtained between the MultiPMomMtch method and the MQS problem at 500 kHz for a source current of 1 A for δ₁ = 0.2, δ₂ = 2.5, and μ₁ = 2000 (left) and for δ₁ = 0.2, δ₂ = 0.1, and μ₁ = 2000 (right).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 12. Distribution of error (|J<sub>PMOR</sub> - J<sub>MQS</sub>|) on current density obtained between the interpolation of projection subspaces method with structure preserving and the MQS problem at 500 kHz for a source current of 1 A for δ₁ = 0.2, δ₂ = 2.5, and μ₁ = 2000 (left) and for δ₁ = 0.2, δ₂ = 0.1, and μ₁ = 2000 (right).</figcaption></figure>

4.5. Hybrid method

We could see in the previous section that the MultiPMomMtch method leads to more accurate results than the interpolation of projection subspaces. However, to apply the MultiPMomMtch method, an affine decomposition is required. If the parameters are material characteristics, the affine decomposition can be obtained easily as shown in the EQS application. This is not so straightforward for geometric parameters. Indeed, for both methods, to manage reduced bases of the same size, the number of degrees of freedom of the mesh should be kept constant when varying the parameters. To reduce numerical noise, the mesh should only be deformed without connectivity modification between nodes and elements. For the multi-moment matching method, there is an additional constraint. Indeed, each element influenced by a geometric deformation will have a new Jacobian matrix which would change the local FEM matrix. If the variation of Jacobian matrix is different for each element, the affine decomposition will have as many parameters as the number of elements. If this number is too large, the application of the multi-moment matching will be intractable (see Section 2.1). In the case of the planar inductor, the geometric deformation is not so complex so we could find a way to keep the variation of the Jacobian matrix identical. This is not always an easy task and it is also highly problem dependent. This is why, in this section, we propose a hybrid method taking advantage of both methods that is presented in the following. This approach mitigates the issue related to the determination of the affine decomposition in the case of the MultiPMomMtch method. We apply first the multi-moment matching method for the material parameter considering different geometry configurations, which would generate several snapshot reduced bases on the offline stage. In the online stage, the interpolation of projection subspaces method is used to interpolate the precomputed reduced bases accounting for new geometric parameter. In this way, we can ignore the affine decomposition in function of geometric parameter required by the MultiPMomMtch method. In the following, we will apply and test the hybrid method to the planar inductor.

As mentioned previously, it is easy for MQS (but also EQS) problems to have an affine decomposition for the material parameter considering that the parameters δ₁ and δ₂ are fixed. In the case of the planar inductor, we have:

$$
\mathbf{K}_A(\delta_1, \delta_2, \mu_1) = \mathbf{K}_\delta(\delta_1, \delta_2) + \mu_1 \mathbf{K}_7 \quad (83)
$$

where the matrix function K<sub>δ</sub>(δ<sub>1</sub>, δ<sub>2</sub>) encapsulates the geometry dependence. With (73), if we fix the parameters (δ<sub>1</sub>, δ<sub>2</sub>) the matrix
K<sub>δ</sub>(δ<sub>1</sub>, δ<sub>2</sub>) is fixed as well. Eq. (83) becomes an affine decomposition. We can then apply the MultiPMomMtch method for the
parameter (s, μ<sub>1</sub>), which would generate a reduced base V<sub>1</sub> for the fixed geometry parameter. If we repeat the process K times for
<figure><img src="image_16.png" /><figcaption>Fig. 13. Log of error of the resistance (left) and the inductance (right) for the hybrid method, the pure interpolation of projection subspaces method with structure preserving.</figcaption></figure>

<table><caption>**Table 6**
Dimension of reduced bases obtained by the hybrid method and the interpolation of projection subspaces method.</caption><thead><tr><th></th><th>Size</th><th>Number of bases</th></tr></thead><tbody><tr><td>Base V for the hybrid method</td><td>1 952 932 × 12</td><td>4</td></tr><tr><td>Base W(p) Interpolation of projection subspaces</td><td>1 952 932 × 16</td><td>8</td></tr></tbody></table>

different sets of parameters ($\delta_1, \delta_2$), we obtain K reduced bases. In our case, we consider K = 4 for the following sets of parameters:

$$ 
\begin{aligned}
\delta_1 &\in \{0.1, 1\}, \\
\delta_2 &\in \{0.1, 2.5\}.
\end{aligned} 
\qquad (84) $$

Once the K = 4 bases are generated in the offline stage, apply the interpolation of projection subspaces method to generate a new basis for each new geometry configuration during the online stage. In this way, we do not rely on an affine decomposition for the geometry parameter. For each geometry configuration, l = 4 expansion points are chosen at the vertices of the cube formed by the domain of variation of the parameters $\mu_1$ and s:

$$ 
\begin{aligned}
\mu_1 &\in [200, 2000], \\
s &\in [2\pi \times 10^3, 10\pi \times 10^5].
\end{aligned} 
\qquad (85) $$

When the MultiPMomMtch method is applied for the affine decomposition (83), a reduced basis composed of 12 vectors is obtained. We can see that the size of the full problem solved in the offline stage is smaller than when the MultiPMomMtch method is applied since the number of parameters has been reduced (1 instead of 3). To evaluate the hybrid method and to compare it with the interpolation of projection subspaces method, we computed the K = 8 projection subspaces for the parameter set:

$$ 
\begin{aligned}
\delta_1 &\in \{0.1, 1\}, \\
\delta_2 &\in \{0.1, 2.5\}, \\
\mu_1 &\in \{200, 2000\}.
\end{aligned} 
\qquad (86) $$

Each projection subspace is still constructed using 8 different frequencies solutions similar to the previous section. The dimensions of the bases obtained by the hybrid method and the interpolation of projection subspaces method are shown in Table 6. We can see that the number of MQS problems solved during the offline stage for the two methods is not the same, the set of parameters $\delta_1, \delta_2$ and $\mu_1$ chosen in the offline stage are the same. During the online stage, the interpolation will not be applied on the same kind of reduced basis. For the hybrid method, the reduced basis is constructed from the MultiPMomMtch method considering material parameters. For the interpolation of projection subspaces method, the reduced basis is obtained by applying POD on snapshots spanning the frequency domain for fixed material and geometry parameters.

The comparison of the interpolation of projection subspaces method and the hybrid method for the global quantities is shown in Fig. 13. The distribution of error of resistance obtained by the hybrid method is shifted towards the left meaning that we have
<figure><img src="image_17.png" /><figcaption>**Fig. 14.** Distribution of error ($|J_{PMOR} - J_{MQS}|$) on current density obtained between the hybrid method and the MQS problem at 500 kHz for a source current of 1 A and $\delta_1 = 0.2$, $\delta_2 = 2.5$, and $\mu_1 = 2000$ (left) and for $\delta_1 = 0.2$, $\delta_2 = 0.1$, and $\mu_1 = 2000$ (right).</figcaption></figure>

improved the accuracy. The distribution of the error of the current density between the hybrid method and the MQS problem is shown in Fig. 14. We can see that the hybrid method enables to improve also the accuracy on local quantities.

**5. Conclusion**

In this paper, the MultiPMomMtch method and the interpolation of projection subspaces method have been applied. These methods have been compared in the context of the DAE associated with the EQS and MQS problem on two realistic examples. The material characteristic parameters as well as the geometry parameters have been considered. Both methods can speed up the construction of the reduced problem when parameters change. The MultiPMomMtch method is faster in the online stage and seems to be more accurate. One major constraint of this method is that it requires an affine decomposition of the parametric problem. We have shown in our application that this decomposition can be obtained even with the geometric parameter. However, this decomposition can be difficult to obtain for more complex deformations of the geometry. This is why we have proposed a hybrid approach mixing MultiPMomMtch and the interpolation of projection subspaces methods. To apply this method, there is no need for an affine decomposition as a function of the geometry parameters. It is shown that this hybrid method is more accurate than the interpolation of projection subspaces method.

The future research consists of applying the PMOR method for the nonlinear problem and on more complex examples such as electrical motors. The reduced problem can be used to couple with circuit equations directly. The equivalent circuit of the reduced problem can be also synthesized, which makes the coupling of the external circuit easier in the electronic circuit simulator.

**CRediT authorship contribution statement**

**Wei Chen:** Writing – original draft, Visualization, Methodology, Conceptualization. **Thomas Henneron:** Writing – review & editing, Supervision. **Stéphane Clénet:** Writing – original draft, Methodology, Conceptualization.

**Declaration of competing interest**

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: This work was supported by EE4.0 project.
Acknowledgments

This work has been achieved within the framework of EE4.0 project (Energie Electrique 4.0). EE4.0 is co-financed by European Union with the financial support of the European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-France.

Appendix A. Implicit PMOR for multi-moment matching method

The algorithm 1 of the implicit PMOR for multi-moment matching method is proposed in [8].

**Algorithm 1** Implicit PMOR for multi-moment matching method [8]

1: Input: matrix $\hat{Q}$, integer $p$, matrices $Q_1, \dots, Q_p$, integer $q$, F vector for SISO problems and matrix for MIMO problems.

2: Output: A matrix V whose column vector is denoted as: $v_1, v_2, \dots$.

3: Initialize $a_1 = 0, a_2 = 0, \text{sum} = 0$.

4: Compute $\mathbf{R}_0 = [\hat{Q}]^{-1}\mathbf{F}$.

5: Compute the first column in V: $v_1 = \mathbf{R}_0 / \|\mathbf{R}_0\|$.

6: Compute the orthonormal columns in $\mathbf{R}_1, \dots, \mathbf{R}_q$ iteratively as below

7: **for** $i = 1, 2, \dots, q$ **do**

8: &nbsp;&nbsp;&nbsp;&nbsp; $a_2 = \text{sum}$;

9: &nbsp;&nbsp;&nbsp;&nbsp; **for** $t = 1, 2, \dots, p$ **do**

10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $a_1 = a_2$ **then**

11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; &nbsp;&nbsp;&nbsp;&nbsp; stop

12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **else**

13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $j = a_1 + 1, \dots, a_2$ **do**

14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{w} = [\hat{Q}]^{-1}\mathbf{Q}_j v_j$;

15: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\text{col} = \text{sum} + 1$;

16: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **for** $k = 1, 2, \dots, \text{col} - 1$ **do**

17: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $h = \mathbf{v}_k^T \mathbf{w}$; $\mathbf{w} = \mathbf{w} - h \mathbf{v}_k$;

18: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end for**

19: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **if** $\|\mathbf{w}\|_2 > \epsilon$ **then**

20: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{v}_{\text{col}} = \frac{\mathbf{w}}{\|\mathbf{w}\|_2}$; $\text{sum} = \text{col}$;

21: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end if**

22: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end for**

23: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end if**

24: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end for**

25: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $a_1 = a_2$;

26: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; **end for**

27: Orthogonalize the columns in V by MGS

Appendix B. Affine decomposition of $k_e$ for partition meshed with prism elements

The dependence on the material parameter $\mu_1$ of $\mathbf{K}_A(\delta_1, \delta_2, \mu_1)$ is similar to the EQS problem. The change of parameter $\delta_1$ and $\delta_2$ corresponds to compression and dilation of bloc 2 and bloc 3 (see Fig. 8), respectively. The matrix is assembled for the five blocs:

$$ \mathbf{K}_A(\delta_1, \delta_2, \mu_1) = \sum_{i=1}^{5} \mathbf{K}_{b,i}(\delta_1, \delta_2, \mu_1) \qquad (B.1) $$

where the matrix $\mathbf{K}_{b,i}(\delta_1, \delta_2, \mu_1)$ is assembled at bloc i. Since both bloc 2, bloc 3 and bloc 5 are meshed with prism elements, the matrices $\mathbf{K}_{b,2}(\delta_1, \delta_2, \mu_1)$, $\mathbf{K}_{b,3}(\delta_1, \delta_2, \mu_1)$ and $\mathbf{K}_{b,5}(\delta_1, \delta_2, \mu_1)$ have a similar form. In the following, we present the construction associated with one of these three matrices and the parameter $\delta$ denotes one fifth of (since those blocs are meshed with 5 layers of prism elements) $\delta_1, \delta_2$, or $2.6 - \delta_2$.

The real element can be defined as the image by a geometric transformation $\Phi_K$ (see Fig. B.15) of a reference element for which the interpolation functions are known. The expressions of the elementary matrices show differential operators applied to the interpolation functions. We have the affine transformation:

$$ \begin{gathered} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} a_2 - a_1 & a_3 - a_1 & 0 \\ b_2 - b_1 & b_3 - b_1 & 0 \\ 0 & 0 & \delta/2 \end{bmatrix} \begin{bmatrix} \hat{x}_1 \\ \hat{x}_2 \\ \hat{x}_3 \end{bmatrix} + \begin{bmatrix} a_1 \\ b_1 \\ c_1 + \delta/2 \end{bmatrix} \\ = \mathcal{J} \begin{bmatrix} \hat{x}_1 \\ \hat{x}_2 \\ \hat{x}_3 \end{bmatrix} + \begin{bmatrix} a_1 \\ b_1 \\ c_1 + \delta/2 \end{bmatrix} \end{gathered} \qquad (B.2) $$
<figure><img src="image_18.png" /><figcaption>Fig. B.15. Reference element to true element.</figcaption></figure>

where $J$ is the Jacobian matrix of the transformation. By denoting that $J_s = \begin{bmatrix} a_2 - a_1 & a_3 - a_1 \\ b_2 - b_1 & b_3 - b_1 \end{bmatrix}$ is the Jacobian matrix for a triangular element, the Jacobian matrix $J$ is written as follows:

$$J = \begin{bmatrix} J_s & 0 \\ 0 & \delta/2 \end{bmatrix} \qquad (B.3)$$

Group the edge element in reference coordinate into matrix form:

$$\hat{w}^1 = [\hat{w}_1^1, \hat{w}_2^1, \dots, \hat{w}_9^1] \qquad (B.4)$$

Let $w^1$ denote the matrix of edge elements in the real coordinate system. We have the following relation:

$$\begin{gathered} w^1 = J^{-1} \hat{w}^1 \\ \text{curl} w^1 = \frac{J}{\det(J)} \text{curl} \hat{w}^1 \end{gathered} \qquad (B.5)$$

Then, the elementary matrix is computed as:

$$k^e = \int_D [\text{curl} w^1]^T \text{curl} w^1 dD \qquad (B.6)$$

The integration can be rearranged as:

$$k^e = \int_0^1 \int_0^1 \int_{-1}^1 [\text{curl} \hat{w}^1]^T \frac{J^T J}{(\det(J))^2} \text{curl} w^1 \det(J) d\hat{x}_1 d\hat{x}_2 d\hat{x}_3 \qquad (B.7)$$

By splitting the Jacobian matrix, one can obtain an affine decomposition of $k^e(\delta)$:

$$k_1^e(\delta) = \int_0^1 \int_0^1 \int_{-1}^1 [\text{curl} \hat{w}^1]^T \frac{[J_s^T J_s \quad 0 \quad 0]}{\det(J_s) \delta/2} \text{curl} w^1 d\hat{x}_1 d\hat{x}_2 d\hat{x}_3 \qquad (B.8)$$

$$k_2^e(\delta) = \int_0^1 \int_0^1 \int_{-1}^1 [\text{curl} \hat{w}^1]^T \frac{\begin{bmatrix} 0 & 0 \\ 0 & \delta^2/4 \end{bmatrix}}{\det(J_s) \delta/2} \text{curl} w^1 d\hat{x}_1 d\hat{x}_2 d\hat{x}_3 \qquad (B.9)$$

$$k^e(\delta) = k_1^e(\delta) + k_2^e(\delta) \qquad (B.10)$$

It can be seen that $k_1^e(\delta)$ and $k_2^e(\delta)$ are proportional to the inverse of $\delta$ and to $\delta$, respectively.

The integrals can be numerically integrated by means of Gaussian quadrature. The global matrix is obtained by assembling the elementary matrices. For example, the matrix for the second bloc is:

$$K_{b,2}(\delta_1, \delta_2, \mu_1) = \bigcup_{e=1}^{N_{ele}} k^e(\delta_1) \qquad (B.11)$$

where $N_{ele}$ is the number of elements in the second bloc.
# Data availability

The authors do not have permission to share data.

# References

[1] L. Sirovich, Turbulence and the dynamics of coherent structures. III. Dynamics and scaling, Quart. Appl. Math. 45 (3) (1987) 583-590, http://dx.doi.org/10.1090/qam/910464.

[2] A. Odabasioglu, M. Celik, L. Pileggi, PRIMA: Passive reduced-order interconnect macromodeling algorithm, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 17 (8) (1998) 645-654, http://dx.doi.org/10.1109/43.712097.

[3] P. Benner, S. Gugercin, K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Rev. 57 (4) (2015) 483-531, http://dx.doi.org/10.1137/130932715.

[4] U. Baur, P. Benner, B. Haasdonk, C. Himpe, I. Martini, M. Ohlberger, Chapter 9: Comparison of methods for parametric model order reduction of time-dependent problems, in: Model Reduction and Approximation, pp. 377-407, http://dx.doi.org/10.1137/1.9781611974829.ch9.

[5] P. Benner, S. Grivet-Talocia, A. Quarteroni, G. Rozza, W. Schilders, L.M. Silveira (Eds.), System- and Data-Driven Methods and Algorithms, De Gruyter, 2021, http://dx.doi.org/10.1515/9783110498967.

[6] P. Benner, W. Schilders, S. Grivet-Talocia, A. Quarteroni, G. Rozza, L. Miguel Silveira (Eds.), Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms, De Gruyter, 2020, http://dx.doi.org/10.1515/9783110671490.

[7] P. Benner, W. Schilders, S. Grivet-Talocia, A. Quarteroni, G. Rozza, L. Miguel Silveira (Eds.), Model Order Reduction: Volume 3 Applications, De Gruyter, 2020, http://dx.doi.org/10.1515/9783110499001.

[8] P. Benner, L. Feng, A robust algorithm for parametric model order reduction based on implicit moment matching, in: A. Quarteroni, G. Rozza (Eds.), Reduced Order Methods for Modeling and Computational Reduction, Springer International Publishing, Cham, 2014, pp. 159-185, http://dx.doi.org/10.1007/978-3-319-02090-7-6.

[9] L. Daniel, O.C. Siong, L. Chay, K.H. Lee, J. White, A multiparameter moment-matching model-reduction approach for generating geometrically parameterized interconnect performance models, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 23 (5) (2004) 678-693, http://dx.doi.org/10.1109/TCAD.2004.826583.

[10] L. Feng, E. Rudnyi, J. Korvink, Preserving the film coefficient as a parameter in the compact thermal model for fast electrothermal simulation, IEEE Trans. Comput.-Aided Des. Integr. Circuits Syst. 24 (12) (2005) 1838-1847, http://dx.doi.org/10.1109/TCAD.2005.852660.

[11] P. Benner, T. Damm, Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic systems, SIAM J. Control Optim. 49 (2) (2011) 686-711, http://dx.doi.org/10.1137/09075041X.

[12] D. Amsallem, C. Farhat, Interpolation method for adapting reduced-order models and application to aeroelasticity, AIAA J. 46 (7) (2008) 1803-1813, http://dx.doi.org/10.2514/1.35374.

[13] D. Amsallem, C. Farhat, An online method for interpolating linear parametric reduced-order models, SIAM J. Sci. Comput. 33 (5) (2011) 2169-2198, http://dx.doi.org/10.1137/100813051.

[14] N.T. Son, T. Stykel, Model order reduction of parameterized circuit equations based on interpolation, Adv. Comput. Math. 41 (5) (2015) 1321-1342, http://dx.doi.org/10.1007/s10444-015-9418-z.

[15] U. Baur, P. Benner, Model reduction for parametric systems using balanced truncation and interpolation, At-Automatisierungstechnik 57 (8) (2009) 411-419.

[16] U. Baur, P. Benner, A. Greiner, J. Korvink, J. Lienemann, C. Moosmann, Parameter preserving model order reduction for MEMS applications, Math. Comput. Model. Dyn. Syst. 17 (4) (2011) 297-317, http://dx.doi.org/10.1080/13873954.2011.547658.

[17] R. Torchio, L. Codecasa, L. Di Rienzo, F. Moro, Fast uncertainty quantification in low frequency electromagnetic problems by an integral equation method based on hierarchical matrix compression, IEEE Access 7 (2019) 163919-163932, http://dx.doi.org/10.1109/ACCESS.2019.2952723.

[18] L. Lombardi, Y. Tao, B. Nouri, F. Ferranti, G. Antonini, M.S. Nakhla, Parameterized model order reduction of delayed PEEC circuits, IEEE Trans. Electromagn. Compat. 62 (3) (2020) 859-869, http://dx.doi.org/10.1109/TEMCC.2019.2919909.

[19] M.A. Khattak, D. Romano, G. Antonini, F. Ferranti, Efficient frequency and time-domain simulations of delayed PEEC models with proper orthogonal decomposition techniques, IEEE Access 12 (2024) 168-178, http://dx.doi.org/10.1109/ACCESS.2023.3347193.

[20] R. Dyczij-Edlinger, O. Farle, Algorithmic enhancements, model-order reduction, and multigrid aspects in contemporary finite element implementations, IEEE Trans. Magn. 45 (3) (2009) 1270-1275, http://dx.doi.org/10.1109/TMAG.2009.2012592.

[21] R. Dyczij-Edlinger, O. Farle, Finite element analysis of linear boundary value problems with geometrical parameters, COMPEL - Int. J. Comput. Math. in Electr. Electron. Eng. 28 (4) (2009) 779-794, http://dx.doi.org/10.1108/03321640910958919.

[22] B. Fröhlich, J. Gade, F. Geiger, M. Bischoff, P. Eberhard, Geometric element parameterization and parametric model order reduction in finite element based shape optimization, Comput. Mech. 63 (5) (2019) 853-868, http://dx.doi.org/10.1007/s00466-018-1626-1.

[23] Y. Paquay, O. Brüls, C. Geuzaine, Nonlinear interpolation on manifold of reduced-order models in magnetodynamic problems, IEEE Trans. Magn. 52 (3) (2016) 1-4, http://dx.doi.org/10.1109/TMAG.2015.2477169.

[24] P. Gunupudi, R. Khazaka, M. Nakhla, T. Smy, D. Celo, Passive parameterized time-domain macromodels for high-speed transmission-line networks, IEEE Trans. Microw. Theory Tech. 51 (12) (2003) 2347-2354, http://dx.doi.org/10.1109/TMTT.2003.820169.

[25] Z. Ren, Influence of the RHS on the convergence behaviour of the curl-curl equation, IEEE Trans. Magn. 32 (3) (1996) 655-658.

[26] L. Feng, P. Benner, A new error estimator for reduced-order modeling of linear parametric systems, IEEE Trans. Microw. Theory Tech. 67 (12) (2019) 4848-4859, http://dx.doi.org/10.1109/TMTT.2019.2948858.

[27] P.-A. Absil, R. Mahony, R. Sepulchre, Riemannian geometry of Grassmann manifolds with a view on algorithmic computation, Acta Appl. Math. 80 (2) (2004) 199-220, http://dx.doi.org/10.1023/B:ACAP.0000013855.14971.91.

[28] B.F. Doolin, C.F. Martin, Introduction to Differential Geometry for Engineers, Courier Corporation, 2013.

[29] J. Ferrer, M. García, F. Puerta, Differentiable families of subspaces, Linear Algebra Appl. 199 (1994) 229-252, http://dx.doi.org/10.1016/0024-3795(94) 90351-4.

[30] R. Mosquera, A. Hamdouni, A.E. Hamidi, C. Allery, POD basis interpolation via inverse distance weighting on Grassmann manifolds, Discret. Contin. Dyn. Syst. - S 12 (6) (2019) 1743-1759, http://dx.doi.org/10.3934/dcds.2019115.

[31] W. Chen, T. Henneron, S. Clénet, T. Delagnes, J. Zou, Model order reduction of an electro-quasistatic problem using CLN method, Finite Elem. Anal. Des. 238 (2024) 104185, http://dx.doi.org/10.1016/j.finel.2024.104185.

[32] L. Taylor, X. Margueron, Y. Le Menach, P. Le Moigne, Numerical modelling of PCB planar inductors: Impact of 3D modelling on high-frequency copper loss evaluation, IET Power Electron. 10 (14) (2017) 1966-1974, http://dx.doi.org/10.1049/iet-pel.2017.0086.

[33] H. Zhang, M. Zhang, W. Yuan, An efficient 3D finite element method model based on the T-a formulation for superconducting coated conductors, Supercond. Sci. Technol. 30 (2) (2016) 024005, http://dx.doi.org/10.1088/1361-6668/30/2/024005.

[34] Y. Le Menach, S. Clenet, F. Piriou, Numerical model to discretize source fields in the 3D finite element method, IEEE Trans. Magn. 36 (4) (2000) 676-679, http://dx.doi.org/10.1109/20.877539.