Author: Philipp Jakobs
Introduction
For applications requiring repeated evaluations, such as optimization or uncertainty quantification, the computational cost of assessing the full frequency response of large-scale finite element models becomes prohibitive. Despite their vast capabilities for modeling complex dynamic engineering problems, practitioners often resort to reduced-order modeling or surrogates to address this challenge. Recently, a data-driven rational approximation approach, known as the AAA algorithm, has been proposed to represent a system's frequency response across the frequency domain. Expanding on this, a recent publication[1] has extended the framework to incorporate system parameters as inputs into the rational approximation, termed as the p-AAA approach.
Theory
For both methods the idea is to approximate the transfer function by a rational function. To achieve this, we generate samples at discrete frequency points, as well as the varying parameter for pAAA. The transfer function is then evaluated at the sampled points/pairs. Those points are then used to build the rational approximation by a mixture of approximation and interpolation. The detailed workflow for both methods is described in the following:
AAA-Algorithm
Algorithm 1: AAA (adapted from [1])
1: Given \{s_i\} \text{ and }\{h_{i}\} = \{H(s_i)\}
2: Initialize: k=0
3: Define \tilde{H} = average(h_{i})and set error \leftarrow \frac{||[h_{i}]-[\tilde{H}]||_\infty}{||[h_{i}||_\infty}
4: while error > desired tolerance do
5: Select (s_{\hat{i}}) by the greedy search
6: Update data partitioning
7: k \leftarrow k +1
8: \sigma_k \leftarrow s_i
9: Build Loewner matrix \mathbb{L}_2
10: Solve \min ||\mathbb{L}_2\bold{a}||_2 \text{ s.t. } || \bold{a}||_2 =1
11: Use \bold{a} to update the rational approximant \tilde{H}(s)
12: error \leftarrow \frac{||[h_{i}]-[\tilde{H}(s_i)]||_\infty}{||[h_{i}||_\infty}
13: end while
14: return \tilde{H}
pAAA-Algorithm
Algorithm 2: pAAA (adapted from [1])
1: Given \{s_i\}, \{p_j\}, \text{ and }\{h_{ij}\} = \{H(s_i,p_j)\}
2: Initialize: k=0 \text{ and } q=0
3: Define \tilde{H} = average(h_{ij})and set error \leftarrow \frac{||[h_{ij}]-[\tilde{H}]||_\infty}{||[h_{ij}||_\infty}
4: while error > desired tolerance do
5: Select (s_{\hat{i}},p_{\hat{j}}) by the greedy search
6: Update data partitioning
7: if s_i was not selected in previous iteration then
8: k \leftarrow k +1
9: \sigma_k \leftarrow s_i
10: end if
11: if p_jwas not selected in previous iteration then
12: q \leftarrow q +1
13: \pi_q\leftarrow p_j
14: end if
15: Build Loewner matrix \mathbb{L}_2
16: Solve \min ||\mathbb{L}_2\bold{a}||_2 \text{ s.t. } || \bold{a}||_2 =1
17: Use \bold{a} to update the rational approximant \tilde{H}(s,p)
18: error \leftarrow \frac{||[h_{ij}]-[\tilde{H}(s_i,p_j)]||_\infty}{||[h_{ij}||_\infty}
19: end while
20: return \tilde{H}
Greedy search
Let \tilde{H} be the rational approximant. The next sampling point is determined by finding the current error maximum.
| \sigma_{k+1} = \arg \max_{i=1,...,N-k} \mid H(\hat{\sigma}_i)-\tilde{H}(\hat{\sigma}_i) \mid |
Greedy search
Let \tilde{H} be the rational approximant. The next sampling points \sigma_k \leftarrow s_i and \pi_q\leftarrow p_j are determined by finding the current error maximum.
| (\sigma_{k+1}, \pi_{q+1}) = \arg \max_{i,j} \mid H(\hat{\sigma}_i,\hat{\pi}_j)-\tilde{H}(\hat{\sigma}_i,\hat{\pi}_j) \mid |
Data partitioning
We partition the data
| \{s_1,...,s_N\} = \{\sigma_1,...,\sigma_k\} \cup \{\hat{\sigma}_1,...,\hat{\sigma}_{N-k}\} \\ |
where the set denoted without a hat are be used for interpolation and the set with a hat is approximated.
The corresponding values H(\sigma_i) are denoted with h_i.
Data partitioning
We partition the data
| \{s_1,...,s_N\} = \{\sigma_1,...,\sigma_k\} \cup \{\hat{\sigma}_1,...,\hat{\sigma}_{N-k}\} \\ \{p_1,...,p_M\} = \{\pi_1,...,\pi_q} \cup \{\hat{\pi}_1,...,\hat{\pi}_{M-q}\}\\ |
where the sets denoted without a hat are used for interpolation and the sets with a hat are approximated.
The corresponding values H(\sigma_i,\pi_j) are denoted with h_{ij}.
Loewner Matrix
In the non-parametric case the Loewner Matrix is calculated by
| \left[\begin{matrix}\frac{H(\hat{\sigma}_1)-h_1}{\hat{\sigma}_1-\sigma_1} & \frac{H(\hat{\sigma}_1)-h_2}{\hat{\sigma}_1-\sigma_2}&...& \frac{H(\hat{\sigma}_1)-h_k}{\hat{\sigma}_1-\sigma_k}\\ \frac{H(\hat{\sigma}_2)-h_1}{\hat{\sigma}_2-\sigma_1} & \ddots && \vdots\\ \vdots & \ &\ddots & \vdots\\ \frac{H(\hat{\sigma}_{N-k})-h_1}{\hat{\sigma}_{N-k}-\sigma_1} & \frac{H(\hat{\sigma}_{N-k})-h_2}{\hat{\sigma}_{N-k}-\sigma_2}&...& \frac{H(\hat{\sigma}_{N-k})-h_k}{\hat{\sigma}_{N-k}-\sigma_k}\end{matrix}\right] |
Loewner Matrix
In the parametric case the Loewner matrix is constructed by concatenating three submatrices:
| \mathbb{L}_2 = [\mathbb{L}_{\sigma\hat{\pi}}^T \mathbb{L}_{\hat{\sigma}\pi}^T \mathbb{L}_{\hat{\sigma}\hat{\pi}}^T]^T \in \mathbb{C}^{(MN-kq)\times kq |
Each submatrix is constructed similar to the AAA-case.
The submatrix \mathbb{L}_{\hat{\sigma}\hat{\pi}} contains all the pairs (\hat{\sigma}_\hat{i},\hat{\pi}_\hat{j}) , \mathbb{L}_{\sigma\hat{\pi}} contains all the pairs (\sigma_i,\hat{\pi}_\hat{j}) and \mathbb{L}_{\hat{\sigma}\pi} contains all the pairs (\hat{\sigma}_\hat{i},\pi_j).The reader is refered to 1 for the individual construction.
Rational Approximation
With \beta_{i}=h_{i}\alpha_{i}, the interpolation for the points \sigma_i is ensured, The final rational approximation is then obtained by
| \tilde{H}(s) = \frac{n(s)}{d(s)} = \sum_{i=1}^k \frac{\beta_i}{s-\sigma_i}/\sum_{i=1}^k \frac{\alpha_i}{s-\sigma_i} |
Rational Approximation
With \beta_{ij}=h_{ij}\alpha_{ij}, the interpolation for the pairs \sigma_i is ensured. The final rational approximation is then obtained by
| \tilde{H}(s,p)=\frac{n(s,p)}{d(s,p)}=\sum_{i=1}^k\sum_{j=1}^q\frac{\beta_{ij}}{(s-\sigma_i)(p-\pi_j)}/\sum_{i=1}^k \sum_{j=1}^q \frac{\alpha{ij}}{(s-\sigma_i)(p-\pi_j)} |
Application and Results
The AAA and the pAAA algorithm were applied to a simple cantilever beam example and a more complex geometric structure, the Kelvin Cell. For the varying parameter, the length of the beam and one side length of the cell were chosen.
Timoshenko-Beam
AAA-Results
For the Timoshenko-Beam model the behaviour of the AAA can be made visible. In the following examples 30 samples in the frequency space were taken.
As the algorithm progresses we see how the number of samples that are interpolated increases. After 16 iterations the algorithm has converged. We can express the accuracy as the H2-error: \epsilon_{H_{2\ }}=\int_{-\infty}^{\infty}{\left|\frac{H\left(\omega\right)-\hat{H}\left(\omega\right)}{H\left(\omega\right)}\right|^2d\omega} = 1.8⋅10^{−5}
pAAA - Results
For pAAA the length of the beam is no longer constant, but is used as a varying parameter.
The above Figure shows the H2 error plotted over the varying beam length. For the blue curve with 11 samples in parameter space and 50 samples in frequency space we can see that the error is very small at the points where parameter samples are lying. In between however it goes up to the magnitude of 10^{−1}. Increasing the samples to 21, we can see that the overall approximation gets much bettter, also in between sample points.
Kelvin-Cell
AAA-Results
The plots show how the accuracy may be affected by the number of samples. For 10 samples, the behaviour for the frequency range from 220Hz to 500Hz is captured quite well. However, as there are not enough samples in the lower frequency range, the rational approximant cannot reconstruct the behaviour accurately. Increasing the parameter samples to 20 we get a very accurate result with an error of \epsilon_{H_{2\ }}=\int_{-\infty}^{\infty}{\left|\frac{H\left(\omega\right)-\hat{H}\left(\omega\right)}{H\left(\omega\right)}\right|^2d\omega} =7.1⋅10^{−6}.
| Full-Order Model | Rational approximation | |
|---|---|---|
| FEM-model | 14 sec | 14 sec |
| Evaluate transfer function at samples | - | 1.36 sec |
| Build approximation | - | 0.04 |
| Frequency sweep | 18 sec | 0.004 sec |
| Total | approx. 32 sec | approx 15.4 sec |
pAAA-Results
For the parametric case the results are quite similar as for the Timoshenko beam example. Whilst for p=11 and s=50 the error is still quite high (around 10^{−1}), we can reduce the error drastically by just increasing the frequency samples.
For this set of samples the approximate runtime can be seen in the following Table:
| Full-Order Model | Rational approximation | |
|---|---|---|
| FEM-model | - | 165 sec |
| Evaluate transfer function at samples | - | 35.6 sec |
| Build approximation | - | 6.3 sec |
| Offline Phase | - | approx. 207 sec |
| FEM-model | 14 sec | - |
| Frequency sweep | 18 sec | 0.07 sec |
| Online phase | approx. 32 sec | approx 0.07 sec |
The rational Approximation needs an offline phase where for each parameter sample a FEM model is built and then evaluated at the frequency samples. This is extra time that needs to be considered. However, after that it is much less computational effort to perform an evaluation of the transfer function over the desired frequency range. For this specific example, we can see that the computation time will break even after 7 new evaluations. This will be exceeded in a lot of applications like UQ or Optimisation.
Conclusions & Outlook
In conclusion, the data-driven modeling framework presented in this article offers a promising approach for accurately representing complex systems. Both the AAA algorithm and the pAAA algorithm demonstrate high accuracy for unseen parameters when a sufficient number of samples are available. Furthermore, the computational time required for evaluations can be significantly reduced compared to traditional methods. While the extension of the framework to incorporate additional parameters shows great potential, challenges such as the curse of dimensionality associated with full grid sampling and the difficulty in efficiently storing large datasets must be addressed. Despite these challenges, the advancements made in data-driven modeling pave the way for more efficient and accurate analysis of complex systems in various fields.
References
1 Andrea Carracedo Rodriguez, , Linus Balicki, Serkan Gugercin. "The p-AAA algorithm for data driven modeling of parametric dynamical systems." (2022).