## 2 Proposed Order Reduction Methodology
### Adaptive Scheme  
To start, because of the aforementioned reasons explained in 2.3, to choose the new expansion point, the logarithm midpoint is obtained from the two old expansion points as  
\[f_{\text{mid}}^{s}=10^{(\log(f_{\text{min}}^{s})+\log(f_{\text{max}}^{s}))/2} \tag{8}\]  
where, \(f_{\text{mid}}^{s}\) denotes the midpoint frequency in logarithmic scale, with \(f_{\text{min}}^{s}\) and \(f_{\text{max}}^{s}\) denoting the minimum and maximum frequencies of the two old points, respectively. Second, the \"local error\" is defined as the error between the reduced-order model (ROM) and the full model at frequency, with the \"global error\" as the error between the ROM and the full model over the frequency interval \(\left[f_{\text{min}},\;f_{\text{max}}\right]\). The tolerances of the \"local error\" and \"global errors\" are initially set as \(tol^{l}\) and \(tol^{s}\), respectively. Moreover, one defines a posteriori-error estimator, as given in (9), to determine whether the current reduced model satisfies the accuracy requirements over a certain frequency band or only at one frequency point. Moreover,  
\[E=\max_{j}\frac{\left\|\mathbf{T}(t)-\mathbf{T}_{j}^{j}(t)\right\|_{2}}{\left\|\mathbf{T} (t)\right\|_{2}} \tag{9}\]  
where, \(\mathbf{T}(t)\) and \(\mathbf{T}_{j}^{j}(t)\) denote the time-domain responses of the \(j\)th variable of full- and reduced-order models under an excitation signal, respectively. Parameter \(N_{\max}^{\rm moments}\) is defined as the limit of the number of moments calculated in one expansion point, and parameter \(N_{\max}^{\rm s}\) as the limit of the number of total expansion points used in the MOR procedure. The larger the number of moments at one frequency, the smaller the error between the ROM and the full model at this frequency is. Therefore, once the error of a frequency is smaller than the predefined tolerance, the algorithm stops moment matching at this frequency. In other words, \(N_{\max}^{\rm moments}\) is equal to the smallest number of moments under which the error between the ROM and the full model is smaller than the tolerance at the frequency. Moreover, if the actual number of the expansion points used is larger than \(N_{\max}^{\rm s}\), the tolerance of the \"local error\" (_tol_) is reduced. To facilitate the implementation of the proposed scheme, it is explained step by step as follows. An expansion point is first chosen (here we choose \(s_{1}=f_{\min}\)). The moment-matching procedure at \(s_{1}\) is then activated sequentially until the local error between the ROM and the full model at frequency \(f_{\min}\) is smaller than the predefined tolerance. The global error (\(E\) over the frequency interval \([f_{\min},\ f_{\max}]\)) checking procedure is then followed. If the error is \(E<tol^{\rm s}\), the algorithm stops and outputs the reduced model; otherwise \(f_{\max}\) is selected as the second expansion point. The moment-matching producer is repeated for \(f_{\max}\) until the \"local error\" between the ROM and the high fidelity model is smaller than the tolerance, and the subsequent block moments are aggregated to the basis formed by the previous ones, with an orthogonal normalization then performed to compress the projection basis. The global error \(E\), between the ROM and high-fidelity model over the frequency interval \([f_{\min},\ f_{\max}]\), is then checked. If the error is \(E<tol\), the algorithm stops and outputs the reduced-model; otherwise a middle frequency point \(f_{\rm mid}\) is determined by (8), between \(f_{\min}\) and \(f_{\max}\) selected as the third expansion point. For convenience of explanation, the serial numbers of 1, 3, 2 are used to represent \(f_{\min}\), \(f_{\rm mid}\), and \(f_{\max}\), respectively.