A Modified AAA Algorithm for Learning Stable Reduced-Order Models from Data

Abstract

In recent years, the Adaptive Antoulas-Anderson (AAA) algorithm has been applied extensively for generating Reduced Order Models (ROMs) of large scale Linear Time-Invariant systems starting from measurements of their transfer functions. When used for Model Order Reduction (MOR) of an asymptotically stable system, the ROMs generated applying AAA are not guaranteed to preserve this fundamental property, thus rendering them impractical in many relevant applications. To overcome this issue, we propose a novel algebraic characterization for the stability of ROMs represented under the AAA barycentric structure. We then translate such characterization into a set of convex semidefinite constraints that can be embedded into the AAA optimization routine to explicitly maximize the model accuracy while ensuring its stability. Finally, we generalize the resulting modeling framework to allow for efficient stable MOR of Multi-Input-Multi-Output systems. An extensive set of numerical experiments provides practical evidence for the effectiveness of the proposed approach and compares its performance with those of available state-of-the-art methods.

1 Introduction

Model Order Reduction (MOR) of large-scale dynamical systems [1] is dire for achieving various goals such as employing efficient simulation or devising robust, automatic control laws. Computing Reduced Order Models (ROMs) directly from input–output data is the only viable approach in many application scenarios. When dealing with Linear Time-Invariant (LTI) systems, the data-driven reduction process is most commonly performed via rational approximation of the system transfer function. In the last decades, a large number of algorithms have been proposed to tackle the problem, e.g. via interpolation [3, 21, 31, 42, 46] or optimization [10, 12, 36]. More recently, the Adaptive Antoulas-Anderson (AAA) algorithm [43] has been intensively and successfully applied to solve a vast collection of rational approximation problems [5, 8, 23, 25, 26, 29, 57], including MOR of LTI systems [18, 28]. The approach builds on the interpolatory model structure introduced in [2], by iteratively increasing its complexity based on a greedy error-minimization approach. At each iteration, the model coefficients are chosen to minimize a linearized model-data error function. Then, an additional interpolation point is enforced at the location where the highest inaccuracy occurs. The iteration continues until a desired approximation error on the data values is enforced (based on a user-defined tolerance). The aforementioned workflow is known to be numerically well conditioned and inherently overcomes the problem of model order selection, which affects most of the state-of-the-art methods (for the latter, tailored empirical solutions have been proposed in the past see e.g. [34] for automated order selection in the Vector Fitting (VF) iteration). For these reasons, AAA has the potential to become the method of choice for accuracy-controlled rational interpolation/approximation.

Unfortunately, AAA is also affected by some relevant shortcomings, among which ROM stability enforcement stands out: when the starting LTI system to be reduced is known to be asymptotically stable, The ROM returned by the AAA algorithm is not guaranteed to be the same. This flaw has also been reported for the Loewner Framework (LF), and it is, in general, a problem for most interpolation-, moment matching-, and least squares-based MOR methods [3]. Since in any numerical simulation, spurious numerical instabilities are to be avoided, the lack of stability guarantees prevents in the first place the application of such methods in professional MOR tools for engineering, that almost invariably exploit stability preserving approximation algorithms such as VF [22, 24, 40]. Due to the relevance of the problem, a vast amount of methods have been developed in the last two decades to cope with the need to construct reliable reduced models that are endowed with the stability property. Many of these methods are based on constrained optimization or perturbation approaches, often applied a posteriori once stability violations are identified. A non-exhaustive list of approaches includes the following [32, 33, 38, 45, 48, 49, 54]. Tailored approaches have been proposed to enforce stability in LF or in AAA. These approaches typically rely on unstable pole flipping or truncation, a combination of both [6, 19, 23, 39, 55], or search for the closest stable approximation to the originally unstable interpolant [27, 30].

The goal of this work is to propose a variation of the AAA algorithm that enforces the asymptotic stability of the computed model by construction, while explicitly optimizing its accuracy, hence avoiding potentially inaccurate post-processing steps based on pole flipping or truncation. Our derivations start by showing that under mild assumptions, the barycentric AAA model structure ${\hat{H}}(s) = \frac{N(s)}{D(s)}$ is asymptotically stable if and only if the rational denominator function D(s) has only asymptotically stable zeros, i.e. is minimum phase. In the second place, exploiting the results of [9], we prove that under the working conditions such requirement is equivalent to the enforcement of a relaxed Positive Realness (PR) condition on D(s), known as Almost Strict Positive Realness (ASPR). In this sense, our work can be regarded as a generalization of the sufficient conditions for stability of ROMs in the Sanathanan-Koerner framework previously introduced in [14, 58]. The proposed stability condition translates into a set of bilinear algebraic constraints involving the AAA barycentric weights that, whenever verified, guarantee the stability of the ROM as a byproduct. Our strategy amounts at enforcing explicitly these constraints in the optimization stage of the AAA algorithm so that both accuracy and stability are achieved concurrently. Since the resulting optimization problem is non-convex, we propose an ad-hoc convex relaxation scheme allowing us to solve for the model unknowns via semidefinite programming, in an automated and efficient manner. This concludes our theoretical derivations.

From the practical standpoint, we suggest how to incorporate the proposed optimization routine into the classical AAA iteration for scalar approximation. The resulting algorithm, named stabAAA, is finally generalized for rational matrix approximation, in order to allow for MOR of Multi-Input-Multi-Output (MIMO) systems. The effectiveness and the scalability features of the new approach are evaluated numerically on a number of practical MOR benchmarks selected from the mechanical, electrical, and acoustic domains, and compared with competing methods in terms of performance. The results show that the proposed approach retains the automation and robustness properties of stabAAA, and that enables computing stable ROMs of massively MIMO systems in a fraction of the time required by available state-of-the-art methods, even in case the latter fail due to excessive memory requirements.

2 Notation and Problem Statement

In the following, $j=\sqrt{-1}$ and $s=\sigma + j\omega \in \mathbb {C}$ denote respectively the imaginary unit and the Laplace variable. Lowercase and uppercase italic letters will be reserved for scalar numbers or functions of the s variable respectively (e.g. $z=H(s): \mathbb {C}\rightarrow \mathbb {C}$). Given a complex number $z\in \mathbb {C}$, $z^*$ represents its conjugate. The operators $\texttt {Re}(\cdot )$ and $\texttt {Im}(\cdot )$ return respectively the real and the imaginary parts of their arguments. Lowercase and uppercase bold letters will be used for vectors and matrices respectively (e.g. ${\textbf {x}}\in \mathbb {R}^n$ and ${\textbf {X}}\in \mathbb {R}^{n\times n}$). The matrices ${\textbf {X}}^T$ and ${\textbf {X}}^{\star }$ are the transpose and the Hermitian transpose of ${\textbf {X}}$. The matrix ${\textbf {I}}_n$ is the identity matrix of size n. The expression ${\textbf {J}}\succ {\textbf {K}}$ (${\textbf {J}}\succeq {\textbf {K}}$) means that the matrix ${\textbf {J}}-{\textbf {K}}$ is positive (semi)definite. Given a matrix ${\textbf {X}}$, its (q, p)-th element will be denoted as $x^{q,p}$, and given a matrix function ${\textbf {X}}(s): \mathbb {C}\rightarrow \mathbb {C}^{Q\times P}$ its (q, p)-th element will be denoted as $X^{q,p}(s)$.

Let us consider a Multi-Input-Multi-Output LTI dynamical system $\varvec{\varSigma }$ in descriptor form

$$\begin{aligned} \varvec{\varSigma }: {\left\{ \begin{array}{ll} {\textbf {E}}\,\dot{\varvec{\zeta }}(t) & \!\!\!\!= {\textbf {A}}\, \varvec{\zeta }(t) + {\textbf {B}}\, {\textbf {u}}(t) \\ {\textbf {y}}(t) & \!\!\!\!= {\textbf {C}}\, \varvec{\zeta }(t) + {\textbf {D}}\, {\textbf {u}}(t) \end{array}\right. }, \end{aligned}$$

(1)

with $\varvec{\zeta }(t) \in \mathbb {R}^{N}$ as the state variable, ${{\textbf {u}}}(t) \in \mathbb {R}^P$ as the control input, and ${{\textbf {y}}}(t) \in \mathbb {R}^Q$ is the observed output. Consequently, we have ${\textbf {A}},{\textbf {E}}\in \mathbb {C}^{N \times N}$, ${\textbf {B}}\in \mathbb {C}^{N \times {P}}$, ${\textbf {C}}\in \mathbb {C}^{{Q} \times N}$, ${\textbf {D}}\in \mathbb {R}^{{Q\times P}}$. When the matrix ${\textbf {E}}$ is non-singular, system (1) is equivalent to a state space realization for the system $\varvec{\varSigma }$. When ${\textbf {E}}={\textbf {I}}_N$, we will use the notation $({\textbf {A}}, {\textbf {B}},{\textbf {C}},{\textbf {D}})$ to denote the state space realization resulting from (1). The transfer function associated to (1) is given by ${{\textbf {H}}}(s)\in \mathbb {C}^{Q\times P}$

$$\begin{aligned} {{\textbf {H}}}(s) = {\textbf {C}}(s {\textbf {E}}- {\textbf {A}})^{-1} {\textbf {B}}+{\textbf {D}}, \end{aligned}$$

(2)

and is assumed to verify condition ${\textbf {H}}^*(s)={\textbf {H}}(s^*)$, so that the impulse response of $\varvec{\varSigma }$ is real-valued. In particular, the poles of ${\textbf {H}}(s)$ are either real or appear in complex conjugate pairs. We assume that all the poles of ${\textbf {H}}(s)$ have a strictly negative real part so that $\varvec{\varSigma }$ is asymptotically stable.

In our setting, the system $\varvec{\varSigma }$ is known only in terms of measurements of its transfer function, retrieved over a set of discrete sampling points located along the positive imaginary axis:

$$\begin{aligned} j\lambda _v \in \varGamma =\{j\lambda _1,\cdots , j\lambda _V\} \subset j\mathbb {R}^+, \quad {{\textbf {H}}}_v={{\textbf {H}}(j\lambda _v)}. \end{aligned}$$

(3)

Starting from this data, our objective is to build an asymptotically stable reduced-order LTI system $\hat{\varvec{\varSigma }}$ of (unknown) size $n\ll N$ with transfer function $\hat{{\textbf {H}}}(s)$ fulfilling the following conditions

$$\begin{aligned} {\hat{{\textbf {H}}}^*(s)=\hat{{\textbf {H}}}(s^*)}, \quad \delta (\{\hat{{\textbf {H}}}(j\lambda _v)\},\{{\textbf {H}}_v\})<\epsilon , \end{aligned}$$

(4)

where $\delta$ is a user defined model-data error function and $\epsilon >0$ an assigned error tolerance. More details will be given in Sects. 3.2 and 5.2.

Due to the stability requirements on $\hat{\varvec{\varSigma }}$, all the poles of ${\hat{{\textbf {H}}}(s)}$ must have a strictly negative real part. Therefore, our problem is to design a rational approximation strategy which concurrently

• 1. guarantees a user-prescribed error bound,

• 2. needs no information on the order n of the resulting model function,

• 3. returns a real-valued and asymptotically stable model function ${\hat{{\textbf {H}}}(s)}$.

3 Background: AAA Algorithm for Scalar LTI MOR

The AAA algorithm, as originally introduced in [43], is a greedy scalar rational approximation scheme that iteratively increases the complexity of the rational approximant in order to meet a user-prescribed error tolerance over the available data samples. The main steps of the algorithm, which is intimately connected to the seminal LF [42], are:

• 1. Express the fitted rational approximants in a barycentric representation, which is a numerically convenient way of representing rational functions [11].

• 2. Select the next interpolation (support) points via a greedy scheme by enforcing interpolation at the point where the error at the previous step is maximal.

• 3. Compute the other variables (the so-called barycentric weights) in order to enforce a linearized least squares approximation on the non-interpolated data.

3.1 Model Structure and Realization

The scalar rational approximant $\hat{H}(s) \in \mathbb {C}$ computed by the original AAA algorithm in [43] is based, at any given step of the iteration, upon the barycentric form:

$$\begin{aligned} \hat{H}(s)= \frac{\sum _{i=1}^k \displaystyle \frac{w_i h_i}{s-z_i}}{\sum _{i=1}^k \displaystyle \frac{w_i}{s-z_i}}=\frac{N(s)}{D(s)}, \end{aligned}$$

(5)

where the $z_i \in \mathbb {C}$ are the barycentric nodes (or support points), assumed to be distinct, the $w_i$’s are the barycentric weights, the $h_i$’s are the barycentric values, and N(s), D(s) are both rational functions. It is to be noted that the following interpolation conditions are satisfied, by construction, at the support points (provided that $w_i \ne 0$), i.e.: $\hat{H}(z_i) = h_i$ for all $1\le i \le k$. Typically, the values of the $h_i$’s are evaluations of the unknown transfer function H(s) to be modeled at the $z_i$ locations, i.e., $h_i:= H(z_i)$. In this case, we say that $\hat{H}(s)$ is a rational interpolant, since it interpolates the original transfer function H(s), at chosen support points $z_i$’s, i.e., $\hat{H}(z_i) = H(z_i)$.

For model order reduction applications, the support points are commonly placed on the imaginary axis, so that $z_i=j\lambda _i \in \varGamma$ as defined in (3). In this setting, for the requirement $\hat{H}^*(s)= \hat{H}(s^*)$ to hold, it is necessary to consider the following modified barycentric expression:

$$\begin{aligned} \hat{H}(s)&=\frac{\sum _{i=1}^k \left( \frac{h_i w_i}{s -j \lambda _i} + \frac{(h_i w_i)^*}{s +j\lambda _i} \right) }{\sum _{i=1}^k \left( \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}\right) } = \frac{\sum _{i=1}^k N_i(s)}{\sum _{i=1}^k D_i(s)}, \end{aligned}$$

(6)

where $N_i(s) = \frac{h_i w_i}{s -j \lambda _i} + \frac{(h_i w_i)^*}{s +j\lambda _i}$ and $D_i(s) = \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}$. This structure was first suggested in [43, Sec. 8] (but never used explicitly), and enforced in [30] and [55], but also in other subsequent publications.

The following Lemma provides a closed-form expression for realizing the system $\hat{\varvec{\varSigma }}$ in descriptor form when its transfer function is structured as in (6).

Lemma 1

A realization for the transfer function $\hat{H}(s)$ in (6) is given by

$$\begin{aligned} \begin{aligned} \tilde{{\textbf {A}}}&= \varvec{\varLambda }- \tilde{{\textbf {B}}} \tilde{{\textbf {R}}} \in \mathbb {C}^{(2k+1) \times (2k+1)}, \ \ \tilde{{\textbf {B}}} = \begin{bmatrix} w_{1}&w_{1}^{*}&\cdots&w_{k}&w_{k}^{*}&1 \end{bmatrix}^{\hspace{-0.83328pt}\textsf{T}}\in \mathbb {C}^{(2k+1) \times 1}, \\ \tilde{{\textbf {C}}}&= \begin{bmatrix} h_{1}&h_{1}^{*}&\cdots&h_{k}&h_{k}^{*}&0 \end{bmatrix} \in \mathbb {C}^{1 \times (2k+1)}, \ \tilde{{\textbf {E}}} = \operatorname {diag}(1,1,\ldots ,1,0) \in \mathbb {C}^{(2k+1) \times (2k+1)},\\ \varvec{\varLambda }&= \operatorname {diag}(j\lambda _{1},-j\lambda _{1},\ldots ,j\lambda _{k},-j\lambda _{k},1) \in \mathbb {C}^{(k+1) \times (k+1)}, \ \tilde{{\textbf {R}}} = \begin{bmatrix} 1&1&\cdots&1 \end{bmatrix} \in \mathbb {C}^{1 \times (2k+1)}, \end{aligned} \end{aligned}$$

(7)

i.e., $\tilde{H}(s) = \tilde{{\textbf {C}}} (s \tilde{{\textbf {E}}} - \tilde{{\textbf {A}}})^{-1} \tilde{{\textbf {B}}}$ satisfies the equality $\tilde{H}(s) = \hat{H}(s)$ for all $s \in \mathbb {C}$.

Based on the above lemma, a real-valued realization for $\hat{H}(s)$ can be obtained via standard approaches, e.g., by following [6].

3.2 The Real-Valued AAA Algorithm

In the following, we summarize the AAA algorithm when the model structure is constrained to obey equation (6). Assuming the availability of a data set in form (3), the AAA algorithm builds the final model (6) following a greedy iterative procedure, driven by an error minimization criterion. This iteration is aimed at returning a model which satisfies

$$\begin{aligned} \max _{v=1,\hdots , V} | \hat{H}(j\lambda _v)-h_v|\le \epsilon , \quad j\lambda _v\in \varGamma , \end{aligned}$$

(8)

where $\epsilon >0$ is an admissible approximation error tolerance. By defining the iteration index $\ell$, and denoting with $\hat{H}_{\ell }(s)$ the model obtained at the $\ell$-th iteration, the algorithm is initialized by defining for $\ell = 0$

$$\begin{aligned} \hat{H}_0(s)=V^{-1}\sum _{v=1}^V h_v, \qquad \varGamma ^{(0)}=\varGamma . \end{aligned}$$

(9)

The $\ell$-th iteration starts by determining the sample point $j\lambda _\ell$ belonging to the current test set $\varGamma ^{(\ell -1)}$ over which the maximum error value is attained,

$$\begin{aligned} j\lambda _\ell = \arg \max _{j\lambda \in \varGamma ^{(\ell -1)}} |H(j\lambda )-\hat{H}_{\ell -1}(j\lambda )|, \end{aligned}$$

(10)

and updating the definition of the model and of the test set consequently

$$\begin{aligned} \hat{H}_\ell (s):= \frac{\sum _{i=1}^\ell \left( \frac{h_i^{(\ell )} w_i^{(\ell )}}{s -j \lambda _i} + \frac{(h_i^{(\ell )} w_i^{(\ell )})^*}{s +j\lambda _i} \right) }{\sum _{i=1}^\ell \left( \frac{w_i^{(\ell )}}{s -j\lambda _i} + \frac{(w_i^{(\ell )})^*}{s +j\lambda _i}\right) }, \qquad \varGamma ^{(\ell )}:= \varGamma ^{(\ell -1)}\setminus \{j\lambda _\ell \}. \end{aligned}$$

(11)

The complex weights $w_i^{(\ell )}$ represent the current model unknowns, that are found by enforcing the linearized approximation of $\hat{H}_\ell (s)$ against H(s) over the updated test set $\varGamma ^{(\ell )}$, for all $s\in \varGamma ^{(\ell )}$ as

$$\begin{aligned} \begin{aligned}&\displaystyle \sum _{i=1}^\ell \left( \frac{w_i^{(\ell )}}{s -j\lambda _i} + \frac{(w_i^{(\ell )})^*}{s +j\lambda _i}\right) H(s) \approx \sum _{i=1}^\ell \left( \frac{h_i^{(\ell )} w_i^{(\ell )}}{s -j \lambda _i} + \frac{(h_i^{(\ell )} w_i^{(\ell )})^*}{s +j\lambda _i} \right) \\&\quad \Leftrightarrow \ \ \Big ( \sum _{i=1}^\ell \frac{H(s) - h_i^{(\ell )}}{s - j\lambda _i} w_i^{(\ell )}\Big ) + \Big ( \sum _{i=1}^\ell \frac{H(s) - (h_i^{(\ell )})^*}{s + j\lambda _i} (w_i^{(\ell )})^*\Big ) := E^{(\ell )}(s) \approx 0. \end{aligned} \end{aligned}$$

(12)

Collecting the unknowns in the vector ${\textbf {w}}^{(\ell )} = [w_1^{(\ell )},(w_1^{(\ell )})^*,\ldots ,w_\ell ^{(\ell )},(w_\ell ^{(\ell )})^*]^T$, condition (50) is enforced by solving the following homogeneous least-squares problem

$$\begin{aligned} {\textbf {w}}_\text {opt}^{{(\ell )}}= \arg \min _{|| {\textbf {w}}^{(\ell )}||_2=1} ||\mathbb {L}^{(\ell )} {\textbf {w}}^{(\ell )}||^2_2 \qquad \mathbb {L}^{(\ell )} \in \mathbb {C}^{(V-\ell )\times 2\ell } \end{aligned}$$

(13)

where $\mathbb {L}^{(\ell )}$ is the Loewner matrix and the constraint $||{\textbf {w}}^{(\ell )}||_2=1$ is introduced to rule out the trivial zero solution. The minimizer ${\textbf {w}}_\text {opt}^{(\ell )}$ is recognized as the right singular vector of $\mathbb {L}^{(\ell )}$ associated to its least singular value. The algorithm stops whenever the maximum error (computed according to the cost function in (10)) is below the threshold tolerance $\epsilon$.

From a numerical standpoint, computing the minimizer of problem (13) via SVD does not guarantee that the coefficients $w_i^{(\ell )}$ are obtained as exact complex conjugate pairs. Since this condition is necessary to ensure that $\hat{H}_{\ell }(s)$ is real-valued, the optimization problem (13) must be reformulated accordingly. This can be done by changing the optimization variables and rewriting the optimization problem in terms of real quantities. First, note that the approximation condition (50) is equivalent to

$$\begin{aligned} \texttt {Re}(E^{(\ell )}(s))\approx 0, \quad \texttt {Im}(E^{(\ell )}(s))\approx 0, \qquad \forall s \in \varGamma ^{(\ell )}, \end{aligned}$$

(14)

when both are are enforced in a least-squares sense; then, by observing the conditions in (14) are linear in the unknowns $\alpha _i^{(\ell )}=\texttt {Re}(w_i^{(\ell )})$ and $\beta _i^{(\ell )}=\texttt {Im}(w_i^{(\ell )})$, we can define the following vector of unknowns, ${\textbf {x}}^{(\ell )}=[ \alpha _1^{(\ell )},\beta _1^{(\ell )}\cdots , \alpha _\ell ^{(\ell )},\beta _\ell ^{(\ell )} ]^T$, and enforce (14) by solving the following homogeneous least-squares problem

$$\begin{aligned} {\textbf {x}}_\text {opt}^{{(\ell )}}= \arg \min _{|| {\textbf {x}}^{(\ell )}||_2=1} ||\mathbb {L}_{\text {real}}^{(\ell )} {\textbf {x}}^{(\ell )}||_2, \ \ \text {with} \ \ \mathbb {L}_{\text {real}}^{(\ell )} \in \mathbb {R}^{2(V-\ell )\times 2\ell }, \end{aligned}$$

(15)

where now, all the involved quantities are real-valued and the row dimension of the matrix $\mathbb {L}_{\text {real}}^{(\ell )}$ is twice that of $\mathbb {L}^{(\ell )}$ since the approximation conditions on real and imaginary components are considered separately. Algorithm 1 summarizes the above procedure in a pseudocode format.

Algorithm 1

The modified AAA algorithm with enforced realness.

4 Stability Characterization and Enforcement for AAA Models

The AAA algorithm has proven to be a reliable and efficient method for the computation of ROMs under the interpolatory transfer function representation due to its concurrent features of error minimization and order selection. However, by construction, the algorithm is not guaranteed to return a rational approximation with Hurwitz poles, so the resulting ROM $\hat{\varvec{\varSigma }}$ may exhibit unstable dynamics even when the underlying system $\varvec{\varSigma }$ is asymptotically stable. This condition is likely to occur when the algorithm is not required to fit exactly the underlying data, i.e. when the error threshold $\epsilon$ is not extremely aggressive, a common scenario in engineering applications. When the ROM $\hat{\varvec{\varSigma }}$ is unstable, it becomes basically useless for a number of practical applications, e.g., any time-domain simulation.

When dealing with rational functions in barycentric form, such as (5) or (6), a strict relation holds between the poles of $\hat{H}(s)$ and the zeros of the denominator function D(s). As will be more precisely formalized next, the poles of $\hat{H}(s)$ can be found only at locations where D(s) exhibits a zero. Based on this observation, we propose to guarantee the stability of $\hat{H}(s)$ by placing such zeros over the open left-half complex plane while solving for the model unknowns. Inspired by previous works related to the stability enforcement of rational barycentric models in the Sanathanan-Koerner approximation framework [14, 58], we achieve the goal by studying the positive realness properties of the denominator function D(s). We show that this approach provides an algebraic characterization for the stability of $\hat{H}(s)$ in terms of the barycentric weights $w_i$. We exploit such characterization to formulate a set of constraints on the model unknowns that, under non-restrictive conditions, are necessary and sufficient for $\hat{H}(s)$ being asymptotically stable.

4.1 Stability Characterization

We start by recalling some basic definitions.

Definition 1

(Relative Degree) Let

$$\begin{aligned} F(s)=\frac{f_1(s)}{f_2(s)} \end{aligned}$$

(16)

be a scalar rational function of s and let $f_1$ and $f_2$ be coprime polynomials with $deg ~f_1(s)=d_1$ and $deg ~f_2(s)=d_2$. Then the relative degree of F(s) is $r=d_2-d_1$, while the McMillan degree of F(s) is $m = \max {(d_1,d_2)}$.

Definition 2

(Minimum Phase [37]) Let $f_1(s)$ and $f_2(s)$ in (16) be coprime polynomials. Then F(s) is said to be Minimum Phase (MP) if

$$\begin{aligned} f_1(s)\ne 0 \quad \forall s: \texttt {Re} (s)\ge 0. \end{aligned}$$

(17)

Additionally, a state space system $({\textbf {A}},{\textbf {B}},{\textbf {C}},{\textbf {D}})$ with $F(s)={\textbf {D}}+{\textbf {C}}(s{\textbf {I}}_n-{\textbf {A}})^{-1}{\textbf {B}}$ is said to be MP if it is stabilizable, detectable [1] and (17) holds.

With these definitions, we can already state the following

Proposition 1

Let

$$\begin{aligned} \hat{H}(s)= \frac{\sum _{i=1}^k \displaystyle \frac{w_i h_i}{s-z_i}}{\sum _{i=1}^k \displaystyle \frac{w_i}{s-z_i}}=\frac{N(s)}{D(s)}=\frac{\frac{p(s)}{l(s)}}{\frac{q(s)}{l(s)}} = \frac{p(s)}{q(s)}, \end{aligned}$$

(18)

with $z_i$ distinct, $w_i\ne 0$, $N(s):=\sum _{i=1}^k \displaystyle \frac{w_i h_i}{s-z_i}$, $D(s):=\sum _{i=1}^k \displaystyle \frac{w_i}{s-z_i}$, and with p(s), q(s), l(s) polynomials in s. Then, the following statements hold true

• Whenever D(s) is Minimum Phase, $\hat{H}(s)$ is asymptotically stable.

• Assuming that the polynomials p(s), q(s) share no roots in the closed right-half complex plane, then the poles of $\hat{H}(s)$ are asymptotically stable if and only if D(s) is Minimum Phase.

Proof

First, notice that when $w_i\ne 0$, the function $D(s)={\frac{q(s)}{l(s)}}$ has exactly k poles, hence l(s) is of degree k and is coprime with q(s). Then, according to Definition 2, when all the roots of q(s) have a strictly negative real part, D(s) is MP. Since the poles of $\hat{H}(s)$ are the roots of q(s) (or a subset of them when zero-poles cancellations occur), MP-ness of D(s) suffices to guarantee that $\hat{H}(s)$ is asymptotically stable. This proves the first statement of the proposition.

When p(s) and q(s) share no roots in the closed right-half complex plane, then $\hat{H}(s)$ is asymptotically stable if and only if all the roots of q(s) have a strictly negative real part. Hence, in this case, the MP-ness of D(s) is equivalent to the stability of $\hat{H}(s)$, and the second statement is proved. $\square$

The following represent classical definitions of Positive Real (PR) and Strictly Positive Real (SPR) scalar rational functions.

Definition 3

(Positive Real (PR) function [15]) A rational function F(s) is positive real if

• 1. F(s) has no poles in $\texttt {Re} (s)>0$

• 2. F(s) is real for all positive real s

• 3. $\texttt {Re} (F(s))\ge 0$ for all $\texttt {Re} (s)>0$

Definition 4

(Strictly Positive Real (SPR) function [15]) A rational function $F(s)\in \mathbb {C}$ that is not identically zero for all s is strictly positive real if $F(s-\tau )$ is PR for some $\tau >0$.

A well-known property of SPR transfer functions [15] is that if F(s) is SPR, so is its inverse 1/F(s). Since the SPR condition requires the regularity of the transfer function over the closed right-half complex plane, the property implies that both the poles and the zeros of an SISO SPR transfer function must be asymptotically stable; therefore any SISO SPR transfer function is also MP.

The following result, which was proven in [51], represents a special case of the well-known Positive Real Lemma for SISO transfer functions and provides the link between the state space representation of an LTI system and the positive realness of its transfer function.

Theorem 1

Consider a SISO state space $({\textbf {A}},{\textbf {B}},{\textbf {C}},{\textbf {D}})$ with $F(s)={\textbf {D}}+{\textbf {C}}(s{\textbf {I}}_n-{\textbf {A}})^{-1}{\textbf {B}}$. Suppose that $\det (s{\textbf {I}}_n-{\textbf {A}})$ has only zeros in the open left half complex plane. Suppose $({\textbf {A}},{\textbf {B}})$ is controllable and let $c>0$, ${\textbf {L}}={\textbf {L}}^T\succ 0$ be an arbitrary real symmetric positive definite matrix. Then a real vector ${\textbf {q}}$ and a real matrix ${\textbf {P}}={\textbf {P}}^T\succ 0$ satisfying

$$\begin{aligned} {\left\{ \begin{array}{ll} {\textbf {P}}{\textbf {A}}+ {\textbf {A}}^T {\textbf {P}}= -{\textbf {q}}{\textbf {q}}^T -c {\textbf {L}}\\ {\textbf {P}}{\textbf {B}}- {\textbf {C}}^T=\sqrt{2{\textbf {D}}}{\textbf {q}}, \end{array}\right. } \end{aligned}$$

(19)

exist if and only if F(s) is SPR and c is sufficiently small.

For state-space systems without feedthrough term, we will make use of the following equivalent statement based on a Linear Matrix Inequality (LMI).

Theorem 2

Let $F(s)={\textbf {C}}(s{\textbf {I}}_n-{\textbf {A}})^{-1}{\textbf {B}}$, with $({\textbf {A}},{\textbf {B}},{\textbf {C}})$ fulfilling the assumptions of Theorem 1. Then the conditions

$$\begin{aligned}&{\textbf {P}}{\textbf {A}}+ {\textbf {A}}^T {\textbf {P}}\prec 0 \end{aligned}$$

(20)

$$\begin{aligned}&{\textbf {P}}{\textbf {B}}={\textbf {C}}^T \end{aligned}$$

(21)

are fulfilled with ${\textbf {P}}={\textbf {P}}^T\succ 0$ if and only if F(s) is SPR.

Proof

This comes directly from Theorem 1 by plugging in ${\textbf {D}}=0$ into the second equality, and also by relaxing the first equality since its right-hand side is an arbitrary negative definite matrix. $\square$

The link between the stability of the AAA model structure (6), MP, and SPR transfer functions is established by the following theorem, proven in [9] and reported here for the SISO case.

Theorem 3

Any strictly proper, minimum phase SISO transfer function given by $F(s)={\textbf {C}}(s {\textbf {I}}_n - {\textbf {A}})^{-1} {\textbf {B}}$ with ${\textbf {C}}{\textbf {B}}>0$ can be made SPR via constant output feedback, so that

$$\begin{aligned} \exists g \in \mathbb {R}: G(s)=\frac{F(s)}{1+gF(s)} = {\textbf {C}}(s {\textbf {I}}_n - {\textbf {A}}+g{\textbf {B}}{\textbf {C}})^{-1}{\textbf {B}}\quad \text {is SPR}. \end{aligned}$$

(22)

The above Theorem states that every strictly proper, SISO MP transfer function with ${\textbf {C}}{\textbf {B}}>0$ can generate an SPR transfer function when closed in feedback with a constant output gain. We highlight that the condition ${\textbf {C}}{\textbf {B}}>0$ is necessary to get SPR-ness of the closed-loop state space $({\textbf {A}}-g{\textbf {B}}{\textbf {C}},{\textbf {B}},{\textbf {C}})$, since in order to fulfill (21), for this system it must hold ${\textbf {C}}{\textbf {B}}={\textbf {B}}^T{\textbf {P}}{\textbf {B}}$ with ${\textbf {P}}$ positive definite. Transfer functions satisfying the requirements of Theorem 3 are referred to in the literature as Almost Strictly Positive Real transfer functions, see e.g. [9].

Theorem 3 can be used to check whether the denominator D(s) in (6) is minimum phase. Thus, it provides a tool to characterize the asymptotic stability of the AAA transfer function $\hat{H}(s)$. The following result comes naturally.

Lemma 2

Let

$$\begin{aligned} \hat{H}(s) =\frac{\sum _{i=1}^k \left( \frac{h_i w_i}{s -j \lambda _i} + \frac{(h_i w_i)^*}{s +j\lambda _i} \right) }{\sum _{i=1}^k \left( \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}\right) } = \frac{N(s)}{D(s)}=\frac{\frac{p(s)}{l(s)}}{\frac{q(s)}{l(s)}} = \frac{p(s)}{q(s)}, \ \ D(s)= {\mathcal {C}}(s {\textbf {I}}_{2k} - {\mathcal {A}})^{-1} {\mathcal {B}}, \end{aligned}$$

(23)

where the realization $({\mathcal {A}}, {\mathcal {B}}, {\mathcal {C}})$ is of size 2k; assume that $w_i\ne 0$ and that $\hat{H}(s)$ has $2k-1$ poles. Then $\hat{H}(s)$ is asymptotically stable if and only if

• 1. ${\mathcal {C}}{\mathcal {B}}>0$ and $\exists g\in \mathbb {R}$ such that

  $$\begin{aligned} G^+(s)= \frac{D(s)}{1+gD(s)} = {\mathcal {C}}(s {\textbf {I}}_{2k} - {\mathcal {A}}+g{\mathcal {B}}{\mathcal {C}})^{-1}{\mathcal {B}}\end{aligned}$$

  (24)

  is SPR or,

• 2. ${\mathcal {C}}{\mathcal {B}}<0$ and $\exists g\in \mathbb {R}$ such that

  $$\begin{aligned} G^-(s)= \frac{-D(s)}{1-gD(s)} = -{\mathcal {C}}(s {\textbf {I}}_{2k} - {\mathcal {A}}-g{\mathcal {B}}{\mathcal {C}})^{-1}{\mathcal {B}}\end{aligned}$$

  (25)

  is SPR.

Proof

With reference to the definition

$$\begin{aligned} D(s)=\sum _{i=1}^k \left( \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}\right) =\frac{q(s)}{l(s)}=\frac{\sum _{m=0}^{2k-1} a_m s^m}{\sum _{n=0}^{2k} b_n s^n}, \end{aligned}$$

(26)

we recall that when $\hat{H}(s)$ has $2k-1$ poles, these poles coincide with the roots of the polynomial q(s). Hence this polynomial must have $2k-1$ roots, implying $a_{2k-1}\ne 0$. As proved in [43], when $w_i\ne 0$, $\hat{H}(s)$ has no poles at $\pm j\lambda _i$; hence q(s) and l(s) are coprime and any realization of D(s) of size 2k is minimal. Noticing that $a_{2k-1}=\sum _{i=1}^k(w_i+w_i^*)$ is the sum of the residues of D(s) and that the quantity ${\mathcal {C}}{\mathcal {B}}$ is realization invariant, we see that the equality $a_{2k-1}={\mathcal {C}}{\mathcal {B}}$ holds for every minimal realization since it holds for the Gilbert realization [1]. Hence, when $\hat{H}(s)$ has $2k-1$ poles, it must hold that ${\mathcal {C}}{\mathcal {B}}\ne 0$. Now, when $\hat{H}(s)$ has $2k-1$ poles, p(s) and q(s) are coprime, and, according to Proposition 1, $\hat{H}(s)$ is asymptotically stable if and only if D(s) is MP.

Consider the case in which ${\mathcal {C}}{\mathcal {B}}>0$; due to Theorem 3, D(s) is MP only if there exists g such that $G^+(s)$ is SPR, so Condition 1 is necessary for the stability of $\hat{H}(s)$. To prove sufficiency, notice that the constant output feedback cannot induce pole-zero cancellations in $G^+(s)$ for $g\in \mathbb {R}$. In fact, such cancellation would occur if both the conditions $G^+(s_z)=0$ and $1+gG^+(s_z)=0$ were verified concurrently, which is impossible. Hence, the zeros of $G^+(s)$ coincide with the zeros of D(s). If $G^+(s)$ is SPR, then these zeros have strictly negative real parts, D(s) is MP, and $\hat{H}(s)$ is asymptotically stable. The same arguments can be applied to prove the statement of the theorem when ${\mathcal {C}}{\mathcal {B}}<0$, by simply considering that D(s) is MP if and only if $-D(s)$ is MP. $\square$

Lemma 2 establishes an explicit algebraic link between the state space realization of the denominator function D(s) and the stability of the AAA model $\hat{H}(s)$. In the following section, we show how to exploit this link in order to constrain the AAA algorithm to generate structurally stable models. Figure 1 provides a graphical illustration of the relations between the various properties discussed above.

Fig. 1

Graphical illustration of the various relations between Minimum-Phase (MP), Almost Strictly Positive Real (ASPR), Positive Real (PR), and Strictly Positive Real (SPR) functions. The intersection between PR and the Not-MP set contains PR functions having zeros on $j\mathbb {R}$

4.2 Exact Stability Enforcement

Consider the denominator function D(s) in (23)

$$\begin{aligned} D(s)=\sum _{i=1}^k \left( \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}\right) = {\mathcal {C}}(s {\textbf {I}}_{2k} - {\mathcal {A}})^{-1} {\mathcal {B}}, \end{aligned}$$

(27)

with $w_i\ne 0$ and $j\lambda _i\in j\mathbb {R}^+$. We define the following minimal state space realization for D(s)

$$\begin{aligned} {\mathcal {A}}=\text {blkdiag}[{\mathcal {A}}_1,\cdots ,{\mathcal {A}}_k],\quad {\mathcal {B}}=[{\mathcal {B}}_1,\cdots , {\mathcal {B}}_k]^T, \quad {\mathcal {C}}=[ \alpha _1,\beta _1,\cdots ,\alpha _k,\beta _k ], \end{aligned}$$

(28)

with

$$\begin{aligned} {\mathcal {A}}_i=\begin{bmatrix} 0& \lambda _i\\ -\lambda _i& 0 \end{bmatrix}, \quad {\mathcal {B}}_i=\begin{bmatrix} 2&0 \end{bmatrix}, \quad w_i=\alpha _i+j\beta _i, \quad \forall i=1,\cdots , k. \end{aligned}$$

(29)

We focus on the optimization problem (15) solved at each iteration of the real-valued AAA algorithm. For readability, we take into account a single iteration $\ell =k\ge 1$ and we drop the iteration index $\ell$, rewriting the problem as

$$\begin{aligned} {\textbf {x}}_\text {opt}= \arg \min _{|| {\textbf {x}}||_2=1} ||\mathbb {L}_{\text {real}} {\textbf {x}}||_2=\arg \min _{\Vert {\mathcal {C}}||_2=1} ||\mathbb {L}_{\text {real}} {\mathcal {C}}^T\Vert ^2_2, \end{aligned}$$

(30)

exploiting the fact that ${\mathcal {C}}=[ \alpha _1,\beta _1,\cdots ,\alpha _k,\beta _k ]={\textbf {x}}^T$. The problem can be equivalently rewritten in terms of ${\mathcal {C}}$ dropping the non-convex norm constraint as

$$\begin{aligned} {\mathcal {C}}_\text {opt}= \arg \min _{{\mathcal {C}}}\Vert \mathbb {L}_{\text {real}} ({\mathcal {C}}^T-{\textbf {x}}_\text {opt})\Vert ^2_2, \end{aligned}$$

(31)

where we assume that the ambiguity on the sign of ${\textbf {x}}_\text {opt}$ is removed by requiring ${\textbf {x}}_\text {opt}^T{\mathcal {B}}>0$. When an asymptotically stable approximant $\hat{H}(s)$ with $2k-1$ poles is required, Lemma 2 states that the AAA model obtained by solving (31) is asymptotically stable if and only if the resulting denominator function D(s) can generate a SPR transfer function under constant output feedback. Therefore, we introduce the required constraints into problem (31), in order to guarantee that any of its feasible solutions will lead to a model with the desired stability property.

First, we observe that the length of the vector defining the cost function of (31) can be reduced by performing the state space transformation

$$\begin{aligned} D(s)\Leftrightarrow (\tilde{{\mathcal {A}}},\tilde{{\mathcal {B}}},\tilde{{\mathcal {C}}})=({\textbf {T}}^{-1}{\mathcal {A}}{\textbf {T}},{\textbf {T}}^{-1}{\mathcal {B}},{\mathcal {C}}{\textbf {T}}), \quad \mathbb {L}_{\text {real}}={\textbf {U}}\varvec{\varSigma }{\textbf {V}}^T, \quad {\textbf {T}}={\textbf {V}}\varvec{\varSigma }, \end{aligned}$$

(32)

where $\mathbb {L}_{\text {real}}={\textbf {U}}\varvec{\varSigma }{\textbf {V}}^T$ is the “thin” SVD of $\mathbb {L}_{\text {real}}$, which is assumed to be full column rank. Notice that this condition is practically verified in any realistic model order reduction process. Using this state space realization, the cost function of (31) is rewritten in terms of $\tilde{{\mathcal {C}}}={\mathcal {C}}{\textbf {T}}$ as

$$\begin{aligned} \Vert \mathbb {L}_{\text {real}} ({\mathcal {C}}^T-{\textbf {x}}_\text {opt})\Vert ^2_2= \Vert \varvec{\varSigma }{\textbf {V}}^T({\mathcal {C}}^T-{\textbf {x}}_\text {opt})\Vert _2^2= \Vert \tilde{{\mathcal {C}}}^T-\bar{{\textbf {x}}}\Vert _2^2, \end{aligned}$$

(33)

where $\bar{{\textbf {x}}}=\varvec{\varSigma }{\textbf {V}}^T {\textbf {x}}_\text {opt}$, leading to the equivalent problem

$$\begin{aligned} \min _{{\mathcal {C}}, \tilde{{\mathcal {C}}}}=\Vert \tilde{{\mathcal {C}}}^T-\bar{{\textbf {x}}}\Vert _2^2 \quad \text {subject to}:\;\tilde{{\mathcal {C}}}={\mathcal {C}}{\textbf {T}}. \end{aligned}$$

Notice that the vector defining the cost function of (34) has size 2k, while the vector entering (31) has size $2(V-k)$. Since the number of data points V is much larger than the number of unknowns 2k in every practical application, the variable transformation substantially reduces the size of the problem.

To enforce model stability, we restrict the feasible set of (34) to the set of vectors $\tilde{{\mathcal {C}}}$ for which Condition (1) of Lemma 2 can be verified, and we consider feasible any vector $\tilde{{\mathcal {C}}}$ such that the function $G^+(s)\Leftrightarrow (\tilde{{\mathcal {A}}}-g\tilde{{\mathcal {B}}} \tilde{{\mathcal {C}}},\tilde{{\mathcal {B}}},\tilde{{\mathcal {C}}})$ fulfills the SPR conditions (20), (21) with suitable $g\in \mathbb {R}$. Thus, according to (22), the constrained optimization problem thus reads

$$\begin{aligned} \min _{{\mathcal {C}},\tilde{{\mathcal {C}}},{\textbf {Q}},g}&\Vert \tilde{{\mathcal {C}}}^T-\bar{{\textbf {x}}}\Vert _2^2\end{aligned}$$

(34)

$$\begin{aligned}&\text {subject to:} \nonumber \\&{\textbf {Q}}={\textbf {Q}}^T\succ 0, \quad {\textbf {Q}}\tilde{{\mathcal {B}}}=\tilde{{\mathcal {C}}}^T, \quad \tilde{{\mathcal {C}}}={\mathcal {C}}{\textbf {T}},\nonumber \\&\tilde{{\mathcal {A}}}^T{\textbf {Q}}+{\textbf {Q}}\tilde{{\mathcal {A}}}-g(\tilde{{\mathcal {B}}} \tilde{{\mathcal {C}}})^T{\textbf {Q}}-g{\textbf {Q}}\tilde{{\mathcal {B}}} \tilde{{\mathcal {C}}}\prec 0. \end{aligned}$$

(35)

As the constraints involved in (34) are necessary and sufficient to verify the stability of model structure (6) when it is assumed to have $2k-1$ poles, the optimal solution $\tilde{{\mathcal {C}}}_\text {opt}^T={\textbf {x}}_\text {opt}$ is feasible whenever the unconstrained AAA solution ${\textbf {x}}_\text {opt}$ leads to an asymptotically stable model. When this is not the case, the stability constraints are active and the optimal solution $\tilde{{\mathcal {C}}}^T_\text {opt}$ will not coincide with $\bar{{\textbf {x}}}$. Eliminating the equality constraints in (34), the problem is equivalently rewritten

$$\begin{aligned} \min _{{\textbf {Q}},g}&\Vert {\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}}\Vert _2^2 \nonumber \\&\text {subject to:} \nonumber \\&{\textbf {Q}}={\textbf {Q}}^T\succ 0, \quad \tilde{{\mathcal {A}}}^T{\textbf {Q}}+{\textbf {Q}}\tilde{{\mathcal {A}}}-2g {\textbf {Q}}\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T {\textbf {Q}}\prec 0 . \end{aligned}$$

(36)

4.3 Stability Enforcement Based on Convex Optimization

Unfortunately, the optimization problem (36) is non-convex, since it requires the fulfillment of a nonlinear matrix inequality in variables ${\textbf {Q}}$ and g. Although a number of methods based on convex relaxations have been proposed to tackle this kind of problem (see e.g. [56]), none of them guarantees the recovery of the global solution. On the other hand, treating the problem as is, employing for example global optimization techniques, would be practically feasible only in case very few optimization variables are involved. In the following, we show that problem (36) can be relaxed to a convex problem that preserves the exactness of the model stability constraints, at the price of a slight modification of the target cost function.

We start by noticing that the non-linear matrix inequality involved in the problem can be turned into a linear one by applying a congruence transformation. Defining ${\textbf {Y}}={\textbf {Y}}^T= {\textbf {Q}}^{-1}$ we have

$$\begin{aligned} \begin{aligned} \tilde{{\mathcal {A}}}^T{\textbf {Q}}+{\textbf {Q}}\tilde{{\mathcal {A}}}-2g {\textbf {Q}}\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T {\textbf {Q}}\prec 0 \quad&\iff \quad {\textbf {Y}}(\tilde{{\mathcal {A}}}^T{\textbf {Q}}+{\textbf {Q}}\tilde{{\mathcal {A}}}-2g {\textbf {Q}}\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T {\textbf {Q}}){\textbf {Y}}\prec 0\\&\iff \quad {\textbf {Y}}\tilde{{\mathcal {A}}}^T+\tilde{{\mathcal {A}}} {\textbf {Y}}-2g\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T\prec 0. \end{aligned} \end{aligned}$$

(37)

Second, recalling that scaling the objective function with an arbitrary positive constant $\eta >0$ does not change the minimizer of the problem, we rewrite (36) equivalently as

$$\begin{aligned} \min _{{\textbf {Q}},{\textbf {Y}},g}&\, \eta \Vert {\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}}\Vert _2^2 \nonumber \\&\text {subject to:} \nonumber \\&{\textbf {Q}}={\textbf {Q}}^T \succ 0, \quad {\textbf {Y}}={\textbf {Q}}^{-1}, \nonumber \\&{\textbf {Y}}\tilde{{\mathcal {A}}}^T+\tilde{{\mathcal {A}}} {\textbf {Y}}-2g\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T\prec 0. \end{aligned}$$

(38)

Although the matrix inequalities that guarantee the model stability are linear in the unknown matrices ${\textbf {Y}}$,${\textbf {Q}}$, the problem is still non-convex due to the presence of the constraint ${\textbf {Y}}={\textbf {Q}}^{-1}$. Applying an inverse Schur complement, the problem is turned into epigraph form as

$$\begin{aligned} \min _{{\textbf {Q}},{\textbf {Y}},g,r}&\,r \end{aligned}$$

(39)

$$\begin{aligned}&\text {subject to:} \nonumber \\&{\textbf {Q}}={\textbf {Q}}^T \succ 0,\quad {\textbf {Y}}={\textbf {Q}}^{-1},\nonumber \\&{\textbf {Y}}\tilde{{\mathcal {A}}}^T+\tilde{{\mathcal {A}}} {\textbf {Y}}-2g\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T\prec 0, \nonumber \\&\begin{bmatrix} r & ({\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}})^T \\ {\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}} & \frac{1}{\eta }{\textbf {I}}_{2k} \end{bmatrix} \succeq 0 . \end{aligned}$$

(40)

The last matrix inequality can be equivalently restated in terms of ${\textbf {Y}}$ by applying the following congruence transformation

$$\begin{aligned} \begin{bmatrix} 1& {\textbf {0}}^T\\ {\textbf {0}}& {\textbf {Y}}\end{bmatrix} \begin{bmatrix} r & ({\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}})^T \\ {\textbf {Q}}\tilde{{\mathcal {B}}}-\bar{{\textbf {x}}} & \frac{1}{\eta }{\textbf {I}}_{2k} \end{bmatrix} \begin{bmatrix} 1& {\textbf {0}}^T\\ {\textbf {0}}& {\textbf {Y}}\end{bmatrix}=\begin{bmatrix} r & (\tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}})^T \\ \tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}} & \frac{1}{\eta }{\textbf {Y}}{\textbf {Y}}\end{bmatrix}\succeq 0. \end{aligned}$$

(41)

Replacing (40) with (41) and introducing the additional constraint ${\textbf {Y}}\succ 0$ allows to eliminate the variable ${\textbf {Q}}$. We then proceed by performing a convex relaxation of (41) relying on the following property, known as Young’s relation [20]

$$\begin{aligned} {\textbf {R}},{\textbf {W}}\in \mathbb {R}^{n\times m}, {\textbf {S}}={\textbf {S}}^T\succ 0 \in \mathbb {R}^{n\times n}\quad \implies \quad {\textbf {R}}^T{\textbf {W}}+ {\textbf {W}}^T{\textbf {R}}\preceq {\textbf {R}}^T {\textbf {S}}^{-1}{\textbf {R}}+{\textbf {W}}^T{\textbf {S}}{\textbf {W}}, \end{aligned}$$

(42)

which, setting ${\textbf {R}}={\textbf {Y}}$, ${\textbf {W}}={\textbf {I}}_{2k}$ and ${\textbf {S}}=\eta {\textbf {I}}_{2k}$ leads to

$$\begin{aligned} \frac{1}{\eta }{\textbf {Y}}{\textbf {Y}}\succeq 2{\textbf {Y}}- {\textbf {I}}_{2k} \eta . \end{aligned}$$

(43)

Due to (43) and (41), any feasible solution of the following problem

$$\begin{aligned} \min _{{\textbf {Y}},g,r}&\,r \nonumber \\&\text {subject to:} \nonumber \\&{\textbf {Y}}\succ 0, \nonumber \\&{\textbf {Y}}\tilde{{\mathcal {A}}}^T+\tilde{{\mathcal {A}}} {\textbf {Y}}-2g\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T\prec 0, \nonumber \\&\begin{bmatrix} r & (\tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}})^T \\ \tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}} & 2{\textbf {Y}}- {\textbf {I}}_{2k} \eta \end{bmatrix} \succeq 0. \end{aligned}$$

(44)

is a feasible solution of (38), and thus of (36), with ${\textbf {Q}}={\textbf {Y}}^{-1}$. Finally, since we assumed that $\eta$ is an arbitrary positive constant, the following implication

$$\begin{aligned} {\textbf {Y}}\succ 0, \quad \begin{bmatrix} r & (\tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}})^T \\ \tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}} & {\textbf {Y}}\end{bmatrix} \succeq 0 \quad \implies \quad \begin{bmatrix} r & (\tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}})^T \\ \tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}}& 2{\textbf {Y}}- {\textbf {I}}_{2k} \eta \end{bmatrix} \succeq 0 \end{aligned}$$

(45)

always holds for sufficiently small $\eta$. Thus, the following relaxation of (36) follows

$$\begin{aligned} \min _{{\textbf {Y}},g,r}&\,r \nonumber \\&\text {subject to:} \nonumber \\&{\textbf {Y}}\succ 0,\nonumber \\&{\textbf {Y}}\tilde{{\mathcal {A}}}^T+\tilde{{\mathcal {A}}} {\textbf {Y}}-2g\tilde{{\mathcal {B}}} \tilde{{\mathcal {B}}}^T\prec 0, \nonumber \\&\begin{bmatrix} r & (\tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}})^T \\ \tilde{{\mathcal {B}}}-{\textbf {Y}}\bar{{\textbf {x}}} & {\textbf {Y}}\end{bmatrix} \succeq 0. \end{aligned}$$

(46)

which is convex in the decision variables and can be solved via standard interior point methods, that are known to be numerically robust and that compute the optimal solution with arbitrary precision in polynomial time [44]. A formal study of the asymptotic complexity of such iterative methods strongly depends on the employed algorithm, implementation, and problem structure. Hence, the time required to solve 46 using a state-of-the-art primal-dual interior point algorithm is experimentally retrieved in Sec. 6.3 for increasing unknowns number.

Once the solution of (46) is computed, the optimal coefficient vector ${\mathcal {C}}$ which defines the AAA model can be recovered exploiting the relations ${\textbf {Y}}={\textbf {Q}}^{-1}$, ${\textbf {Q}}{\tilde{{\mathcal {B}}}}=\tilde{{\mathcal {C}}}^T$, $\tilde{{\mathcal {C}}} = {\mathcal {C}}{\textbf {T}}$. The corresponding solution will guarantee the asymptotic stability of the resulting ROM.

5 Real-Valued AAA Algorithm with Guaranteed Stability

Here, we present a version of the scalar AAA algorithm for model order reduction of scalar LTI systems with stability constraints based on the convex programming formulation (46). Then, we propose our extension of the algorithm for MOR of MIMO systems.

5.1 Embedding Stability Constraints in AAA

Enforcing stability in the classical AAA algorithm could be done by simply modifying Step 4 in Algorithm 1, by solving (46) for the optimal matrix ${\textbf {Y}}_\text {opt}$ in place of (15), and recovering the required coefficient vector ${\textbf {x}}_\text {opt}^{(\ell )}$ via prescribed variable transformations. However, the resulting algorithm at each iteration should solve a semi-definite program in place of simply computing the SVD of $\mathbb {L}_{\text {real}}^{(\ell )}$. Hence, we found that a more efficient and effective approach is to solve (46) only once after the last iteration, and only if the model returned by Algorithm 1 is found to be unstable. In case the resulting stable model does not meet the required maximum error tolerance, we simply perform more AAA iterations and increase the number of poles of $\hat{H}(s)$. A pseudocode for the proposed approach, which we call stable AAA, or stabAAA, is given in Algorithm 2.

Algorithm 2

stabAAA: scalar AAA algorithm with stability enforcement.

5.2 Extension to MIMO Order Reduction

So far, the discussion was focused on deriving a stable version of the scalar AAA algorithm for model order reduction of scalar LTI systems. Here, we show that the proposed approach admits a straightforward extension that allows to perform efficiently model order reduction of multivariable systems. We extend the AAA model structure (6) to the MIMO case, using a representation from [13] (as a modified version of the format mentioned as Bary-A in [29]):

$$\begin{aligned} \hat{{\textbf {H}}}(s)=\frac{\sum _{i=1}^k \left( \frac{{\textbf {H}}_i w_i}{s -j \lambda _i} + \frac{({\textbf {H}}_i w_i)^*}{s +j\lambda _i} \right) }{\sum _{i=1}^k \left( \frac{w_i}{s -j\lambda _i} + \frac{w_i^*}{s +j\lambda _i}\right) } = \frac{{\textbf {N}}(s)}{D(s)} \in \mathbb {C}^{Q\times P},\quad D(s) \in \mathbb {C}, \end{aligned}$$

(47)

where all the symbols retain the same meaning of Sec. 3. The above model structure represents each element of $\hat{{\textbf {H}}}(s)$ according to the scalar (real-valued) model structure representation (6), so that the interpolation condition $\hat{{\textbf {H}}}(j\lambda _i)={\textbf {H}}_i$ holds by construction. Additionally, the barycentric form (47) is defined upon the same barycentric denominator D(s), implying that each element of $\hat{{\textbf {H}}}(s)$ shares the same set of common poles, which coincide with the zeros of D(s). Hence, the stability characterization introduced in Lemma 2 remains valid also for the proposed MIMO model structure (47), since the ROM with transfer function $\hat{{\textbf {H}}}(s)$ is asymptotically stable if and only if all entries of $\hat{{\textbf {H}}}(s)$ have poles with strictly negative real part.

Based on the above considerations, we now introduce a matrix-valued version of Algorithm 2. Being $\ell$ the iteration index, we define the initialization quantities

$$\begin{aligned} \hat{{\textbf {H}}}_0(s)=V^{-1}\sum _{v=1}^V {\textbf {H}}_v, \qquad \varGamma ^{(0)}=\varGamma . \end{aligned}$$

(48)

At any $\ell$-th iteration, the model $\hat{{\textbf {H}}}_\ell (s)$ represents an updated version of $\hat{{\textbf {H}}}_{\ell -1}(s)$:

$$\begin{aligned} \hat{{\textbf {H}}}_\ell (s):= \frac{\sum _{i=1}^\ell \left( \frac{{\textbf {H}}_i w_i^{(\ell )}}{s -j \lambda _i} + \frac{({\textbf {H}}_i w_i^{(\ell )})^*}{s +j\lambda _i} \right) }{\sum _{i=1}^\ell \left( \frac{w_i^{(\ell )}}{s -j\lambda _i} + \frac{(w_i^{(\ell )})^*}{s +j\lambda _i}\right) }, \quad j\lambda _\ell = \arg \max _{j\lambda \in \varGamma ^{(\ell -1)}} \Vert {\textbf {H}}(j\lambda )-\hat{{\textbf {H}}}_{\ell -1}(j\lambda )\Vert _2, \end{aligned}$$

(49)

where the error criterion guiding the selection of $j\lambda _\ell$ can be regarded as a matrix generalization of (10). The test set is then updated as in (11). Defining the linearized matrix error function,

$$\begin{aligned} {\textbf {E}}_\ell (s)= \sum _{i=1}^\ell \left( \frac{w_i^{(\ell )}}{s -j\lambda _i} + \frac{(w_i^{(\ell )})^*}{s +j\lambda _i}\right) {\textbf {H}}(s) - \sum _{i=1}^\ell \left( \frac{{\textbf {H}}_i w_i^{(\ell )}}{s -j \lambda _i} + \frac{({\textbf {H}}_i w_i^{(\ell )})^*}{s +j\lambda _i} \right) , \end{aligned}$$

(50)

the model coefficients $w_i^{(\ell )}$ are found by enforcing in least-squares sense

$$\begin{aligned} {\texttt {Re}}(E_\ell ^{q,p}(s))\approx 0, \; {\texttt {Im}}(E_\ell ^{q,p}(s))\approx 0 \quad \forall s\in \varGamma ^{(\ell )}, \quad p=1,\cdots , P,\,q=1,\cdots , Q, \end{aligned}$$

(51)

being $E_\ell ^{q,p}$ the (q, p)-th element of ${\textbf {E}}_\ell (s)$. Simple inspection shows that each element-wise condition in (51) is equivalent to the scalar approximation condition (14). Hence, the optimal values of the current weights $w_i^{(\ell )}$ are found by solving

$$\begin{aligned} {\textbf {x}}_\text {opt}^{{(\ell )}}= \arg \min _{|| {\textbf {x}}^{(\ell )}||_2=1} ||\mathbb {L}_{M}^{(\ell )} {\textbf {x}}^{(\ell )}||_2, \ \ \text {with} \ \ \mathbb {L}_{M}^{(\ell )} \in \mathbb {R}^{2(V-\ell )PQ\times 2\ell }, \end{aligned}$$

(52)

where the unknown vector ${\textbf {x}}^{(\ell )}$ is defined as in (14), while

$$\begin{aligned} \mathbb {L}_{M}^{(\ell )}=\begin{bmatrix} \mathbb {L}_{\text {real},1,1}^{(\ell )T}{,\,}&\mathbb {L}_{\text {real},1,2}^{(\ell )T}&,\,\cdots&{,\,} \mathbb {L}_{\text {real},q,p}^{(\ell )T},\,\cdots&{,\,} \mathbb {L}_{\text {real},Q,P}^{(\ell )T} \end{bmatrix}^T, \end{aligned}$$

(53)

being $\mathbb {L}_{\text {real},q,p}^{(\ell )}$ the real-valued Loewner matrix computed as for the scalar case (14), with respect to the (q, p)-th element of ${\textbf {E}}_\ell (s)$. The solution to (52) is obtained again via SVD. We remark that solving iteratively problem (52) can be done efficiently, following e.g., the approach described in [13, Sec.V], since the matrix $\mathbb {L}_{M}^{(\ell )}$ is an update of $\mathbb {L}_{M}^{(\ell -1)}$. The iteration stops when the model reaches a prescribed accuracy. As an accuracy measure, in this work we propose to evaluate the Root Mean Square (RMS) error committed by the model against all the available data

$$\begin{aligned} \delta _M^{(\ell )}=\sqrt{\frac{1}{VPQ}\sum _{p=1}^P\sum _{q=1}^Q\sum _{v=1}^V \big (h^{q,p}_v-\hat{H}^{q,p}_{\ell }(j\lambda _v)\big )^2 }, \end{aligned}$$

(54)

where $h^{q,p}_v,\hat{H}^{q,p}_{\ell }(j\lambda _v)$ are the (q, p)-th elements of ${\textbf {H}}_v$ and $\hat{{\textbf {H}}}_\ell (j\lambda _v)$ respectively. Based on (54), the fitting algorithm stops when it holds $\delta _M^{(\ell )}<\epsilon$, being $\epsilon$ the user defined error tolerance. After the final iteration, in case the model is found unstable, stability enforcement is performed as in Algorithm 2: in fact, the approach introduced in Sects. 4.2, 4.3 can be applied seamlessly for the generation of MIMO ROMs, by simply replacing $\mathbb {L}_{\text {real}}$ with $\mathbb {L}_{M}^{(\ell )}$ in (30) and performing the same subsequent derivations. Notice that the size of the semidefinite program (36) does not depend on the number of input or output of the system. Hence, the choice of model structure (47) guarantees that stable MIMO ROMs can be built in a particularly efficient manner, as will be experimentally shown in Sec. 6.

5.3 Can We Always Attain the Desired Accuracy?

We highlight that there are at least two relevant scenarios of practical interest in which a stable model meeting an arbitrary accuracy may not exist. These scenarios occur when

• 1. The target transfer function ${\textbf {H}}(s)$ exhibits one or more unstable poles or, in general, unstable dynamics. Since a closed-form expression for ${\textbf {H}}(s)$ may be unavailable, it is not generally possible to infer from the available data if such a scenario occurs. If it does, a constrained stable rational approximation $\hat{{\textbf {H}}}(s)$ for ${\textbf {H}}(s)$ which satisfies a maximum error bound with arbitrarily small threshold $\epsilon$ may not exist.

• 2. The data samples may include noise and/or some non-causal components, both of which are incompatible with an arbitrarily accurate rational approximation with constrained stability. Such spurious components can come from direct measurement or even from the discretization error of Partial Differential Equations providing a first-principle physical description of the underlying dynamics. Also, models of material properties that are not causal (e.g., as fitted from measurements using non-causal models). These difficulties are well-known and have been extensively studied [54]. Detection of non-causal data components is particularly challenging since any finite set of data samples is compatible with a causal model (e.g. by computing a rational interpolation with order equal to the number of data points). However, severe overfitting occurs when enforcing model causality hence stability when interpolating or approximating non-causal data. Efficient and reliable algorithms exist to detect such situations [52, 53].

In the above scenarios, the approximation error of a constrained stable rational model saturates to a minimum allowed level even when arbitrarily increasing model order [54]. In this work, we assume that the error threshold imposed as a stopping criterion for model order selection is larger than the amount of non-causal components in the data.

6 Numerical Results and Applications

We apply the proposed stability enforcement approach and algorithms to generate ROMs for 6 dynamical systems related to electrical, mechanical, and acoustic domains. We then proceed to compare the results with those obtained by applying other state-of-the-art data-driven MOR approaches. Each dataset is processed upon data standardization, performed as follows:

$$\begin{aligned} f_{Max} = \max _{v} |\lambda _v/2\pi |, \quad h_{Max}= \max _{v,p,q} |h^{q,p}_v |, \quad {\textbf {H}}_v\leftarrow {\textbf {H}}_v/h_{Max}, \quad \lambda _v\leftarrow \lambda _v/f_{Max}. \end{aligned}$$

All the experiments are performed in MATLAB using a workstation equipped with 32 GB of memory and a 3.3 GHz Intel i9-X7900 CPU. To solve problem (46) we use the YALMIP toolbox [41] and the solver MOSEK [4]. A sample implementation for Algorithm 2 can be found at https://github.com/tomBradde/stabAAA.

6.1 Testing the Proposed Approach

6.1.1 A High-Speed Printed Circuit Board Link

We consider a high-speed electrical interconnection between two printed circuit boards (PCB), first presented in [47]. Our objective is to derive a SISO ROM for the output admittance at one port of the interconnect. A high-fidelity characterization for this admittance function is computed via an electromagnetic field solver, as explained in [47]. The available dataset is composed of pairs $(j\lambda _v,h_v)=(j2\pi f_v,h_v), v=1,\cdots , 497, f_v \in [0.06,10]$ GHz. We apply Algorithm 1 with error threshold tolerance $\epsilon =0.001$. The algorithm returns a stable model with 41 poles, meeting the required accuracy after 22 iterations. The transfer function of the model is reported against the available data in the right panel of Fig. 2. The left panel shows a comparison between the stable identified poles and those obtained by removing the stability constraints.

Fig. 2

High Speed PCB link test case. Left panel: zoom on the dominant poles identified by the standard AAA algorithm (hollow black circles) and by stabAAA (red asterisks). Right panel: magnitude plot of the reference data against the response of the model obtained with stabAAA (Color figure online)

6.1.2 International Space Station (ISS) 1R Model

In this section, we analyze a structural model for component 1R (Russian service module) of the ISS. This model was used for MOR purposes in, among others, the work in [35], and has become a standard MOR benchmark in the last two decades. The data for this test case are obtained considering one scalar element of the system transfer matrix, evaluating the full order model (with a state space of dimension $N=270$) along the imaginary axis $(j\lambda _v,h_v), v=1,\cdots , 400, \lambda _v \in [0.1,100]$ rad/s. Applying Algorithm 1 with error threshold tolerance $\epsilon =10^{-4}$ returns a stable model of order 61 that fulfills the required accuracy requirements after 32 iterations. The model response is compared against the data in the right panel of Fig. 3, while the stable poles are compared with those obtained by removing the stability constraint in the left panel.

Fig. 3

ISS test case. Left panel: zoom on the dominant poles identified by the standard AAA algorithm (hollow black circles) and by stabAAA (red asterisks). Right panel: magnitude plot of the reference data against the response of the model obtained with stabAAA (Color figure online)

6.1.3 Acoustic Absorber

The system considered in this section is an acoustic absorber, composed of a cavity with one side covered by a layer of poroelastic material. The structure details are available in [50], and the related data are taken from [7]. The data samples for this example are $(j\lambda _v,h_v)=(j2\pi f_v,h_v), v=1,\cdots , 901, f_v \in [0.1,1]$ kHz. We first apply Algorithm 1 with $\epsilon =10^{-8}$ and $M_{max}=0$. With these settings, an unstable model of order 31 fulfilling the required error is found after $k=16$ iterations. After the stability enforcement is performed, the resulting stable model is affected by an error $\delta =1.92\times 10^{-6}$, that does not meet the accuracy requirements. The results of the fitting and of the stable poles location are show in Fig. 4. As will be shown in Sec. 6.2, similar or worse accuracy is obtained for the competing methods when ROMs with 31 poles are considered.

Fig. 4

Acustic absorber test case. Left panel: zoom on the dominant poles identified by the standard AAA algorithm (hollow black circles) and by stabAAA (red asterisks). Right panel: magnitude plot of the reference data against the response of the model obtained with stabAAA (Color figure online)

As a second experiment, we apply again the proposed algorithm with $M_{max}=5$ and $\theta =0.1$. With these settings, the algorithm stops after $k=17$ iterations, generating a stable model without the need for stability enforcement. The model error metric is $\delta = 3.62\times 10^{-11}$ that meets the prescribed accuracy requirements.

6.1.4 Parallel Metal Plates

The fourth test case is a 2-D distributed system representing a canonical Power Delivery Network (PDN) on a PCB. The system is composed of two parallel metallic plates with length $l_x = 12$ cm and width $l_y = 10$ cm, separated by a dielectric of height $d = 1$ mm. Five ideal electrical ports are defined between the top and the bottom plates, and one port is shunted with a $\mathrm {20~m\varOmega }$ lumped resistor. The inputs of the system are $P=5$ ideal voltage sources placed in correspondence of 5 ports. The outputs are the $Q=5$ related port currents. Hence, the data samples $(j\lambda _v,{\textbf {H}}_v)=(j2\pi f_v,{\textbf {H}}_v), v=1,\cdots , 2000, f_v \in [0.001,2]$ GHz are obtained computing the admittance matrix of the structure via analytical solution of the field equations via modal expansion. See [16, Sec. VI.A] for further details.

We retrieve a ROM for the system applying the MIMO stabAAA algorithm, setting $\epsilon =10^{-5}$. The algorithm stops after 14 iterations returning a ROM transfer matrix with 27 poles. The stability enforcement effectively removes one real unstable pole identified by the unconstrained algorithm counterpart. Fig. 5 reports the error committed by the algorithm over each element of the admittance matrix. In Sec. 6.2 it will be shown how enforcing stability of the AAA model by truncating this unstable pole may have an unwanted effect on the resulting ROM.

Fig. 5

The RMS error committed by the stabAAA model over each scalar element of the sampled transfer matrix for the example of Sec. 6.1.4

6.1.5 A Real Power Delivery Network

This example concerns a real PDN system for a mobile computational platform with $P=Q=144$. This test case was first considered in [17]. The dataset consists of samples of the impedance matrix of the structure $(j\lambda _v,{\textbf {H}}_v)=(j2\pi f_v,{\textbf {H}}_v), v=1,\cdots , 85$, $f_v \in [10^{-8},2.5]$ GHz, retrieved via high fidelity mixed lumped and electromagnetic characterization of the physical structure. Applying the proposed approach with an error tolerance $\epsilon =5\times 10^{-5}$ a stable ROM with the prescribed accuracy is obtained after 13 iterations, and explicit stability enforcement was needed to remove one unstable pole. The left panel of Fig. 6 shows the RMS error of the resulting ROM against each scalar element of the reference data matrix.

Fig. 6

The RMS error committed by the stabAAA models over each scalar element of the sampled transfer matrices for the example of Sec. 6.1.5 (right panel) and Sec. 6.1.6 (left panel)

6.1.6 Large Scale Parallel Metal Plates

The last example is a large-scale MIMO system with $P=Q=289$ inputs and outputs. We consider again the parallel metal plates structure of Sec. 6.1.4 defining in this case with length $l_x = 20$ cm and width $l_y = 10$ cm, a dielectric of height 5 mm, and 289 lumped electrical ports placed in correspondence of a $17\times 17$ grid with nodes uniformly spaced over the plates. We build a ROM for this structure starting from samples of its impedance matrix, considering a dataset $(j\lambda _v,{\textbf {H}}_v)=(j2\pi f_v,{\textbf {H}}_v), v=1,\cdots ,350$, $f_v \in [0.1, 1]$ GHz. With an error threshold $\epsilon =10^{-5}$, the proposed approach returns a stable model after $k=9$ iterations. The stability enforcement removes four unstable poles identified by the unconstrained $\texttt {AAA}.$ The right panel of Fig. 6 reports the RMS error committed by the ROM over each element of the data transfer matrix.

Table 1 Error metrics for the SISO test cases and generated models

Table 2 Error metrics for the MIMO test cases and generated models

6.2 Comparison with State of the Art Methods

The performance of the proposed approach is now compared to state-of-the-art methods in terms of modeling error and runtime. For each of the considered test cases, and for each ROM ${\textbf {X}}(s)$, we compute the following error metrics

$$\begin{aligned} E_{RMS}=\sqrt{\frac{1}{VPQ}\sum _{p=1}^P\sum _{q=1}^Q\sum _{v=1}^V \big (h^{q,p}_v- {\textbf {X}}^{q,p}(j\lambda _v)\big )^2 },\qquad \, E_2=E_{RMS}\sqrt{VPQ} \end{aligned}$$

(55)

and

$$\begin{aligned} E_\infty = \max _{v=1,\ldots ,V} \Vert {\textbf {H}}_v-{\textbf {X}}(j\lambda _v)\Vert _2. \end{aligned}$$

(56)

The competing methods are compared by generating ROMs of the same complexity of that obtained by applying our approach as in Sec. 6.1. The results of this comparison are reported in Tab. 1 (for the SISO case), in Tab. 2 (for the MIMO case), together with the discussion below.

• AAA: the results of this method are obtained via our algorithms without enforcing stability. For all of the considered test cases (except the acoustic absorber case), the proposed approach matches almost exactly the accuracy of the unconstrained counterpart. The time requirement difference is due to the solution to the constrained optimization problem and the YALMIP parser overhead.

• smiAAA: this method was recently introduced in [23] and is the only one that preserves interpolation while enforcing stability via pole flipping. We used the publicly available implementation of the algorithm to perform the experiments. The approach performs similarly or worse (ISS and absorber test-cases) than stabAAA and is significantly less efficient in the case of massively MIMO systems. For the last experiment, results were not reported due to excessive time requirements (longer than 6 h).

• AAA$^2$: this approach was presented in [55] and enforces stability in AAA by truncation of unstable poles. Results show that in some cases this approach can severely degrade the model accuracy (even by orders of magnitude). Also, the interpolation property of the barycentric model structure is lost due to the residue optimization stage that follows the unstable pole truncation.

• rational: this is a MATLAB function for rational approximation. The related documentation refers to [43] so that it is intended that the function is based on AAA. Unfortunately, no further details are provided to the user, and it is unclear how stability is enforced in this algorithm. Also, the ROMs returned by rational do not guarantee interpolation and are less accurate than those obtained via stabAAA except for the ISS test case. For the massively MIMO parallel plates structure, results are not available since the memory requirements of the method exceed our hardware availability.

• VF: we applied an in-house implementation of the VF algorithm for the SISO test cases and the MATLAB function vecfit3() for MIMO examples (available at https://www.sintef.no/en/software/vector-fitting/downloads/vfit3/). The latter implementation is based on the efficient formulation [24]. VF performance is on par with those of the proposed approach, with error mismatches that are irrelevant to any practical extent. From the computational efficiency standpoint, stabAAA grants superior efficiency when processing systems with large input/output counts. Notice that for VF, the reported runtimes do not include any order selection step, which is expected to further increase the time requirements of the method.

• rkfit: we applied the MATLAB implementation of rkfit available at http://guettel.com/rktoolbox/, by requiring stable ROMs with real impulse response. For the sake of performance comparison, the same considerations as for the VF algorithm can be drawn. In Tab. 2 we do not report the results for the last two MIMO applications since the memory requirements of the algorithm exceed the available hardware.

Fig. 7

Error metric and runtime of stabAAA for increasing number of model poles

6.3 stabAAA Performance for Increasing Poles Count

In this last section, we consider again the ISS benchmark, gathering a MIMO dataset $(j\lambda _v,{\textbf {H}}_v)=(j2\pi f_v,{\textbf {H}}_v), v=1,\cdots , 1450$, $f_v \in [0.2, 156]$ Hz. We apply our algorithm with decreasing error tolerance $\epsilon$ within the interval $[0.05,5\times 10^{-8}]$. As a result, we retrieve a collection of ROMs with a number of poles ranging from 3 to 163. The value of $E_{RMS}$ computed with such models is shown in Fig. 7. The figure shows how stabAAA robustly tracks the error committed by the unconstrained AAA also when the pole count exceeds one hundred, confirming the robustness of the proposed optimization framework. In the same figure, the runtime of the algorithm is also reported. The constrained optimization problem is solved in a matter of seconds for a small to moderate number of poles. When the number of poles exceeds one hundred, the time requirement starts to be significant, as could be expected by the polynomial complexity of the employed interior-point optimization solver. At the same time, it is to be noted that requiring ROMs to meet extremely aggressive error tolerances would defeat the purpose and the benefits of the order reduction process per se, which is intended to generate models of small to moderate complexity.

7 Conclusion

In this paper, we proposed a novel formulation of the AAA algorithm which generates certified asymptotically stable ROMs of LTI systems. The main enabling factor to achieve this goal is a set of novel algebraic conditions on the AAA barycentric weights that, whenever verified, guarantee the stability of the model by construction. Since these conditions are non-convex (on the model coefficients), we propose to enforce them during model generation. This is accomplished by solving a relaxed semi-definite programming problem that can be efficiently handled by any off-the-shelf convex optimization solvers. The complexity of such a problem does not depend on the input–output dimension of the underlying system. The resulting algorithm is experimentally proved to be fully automated and robust, and particularly efficient when applied for the generation of ROMs with a large number of inputs and outputs and a small-to-moderate number of poles. As for similar methods based on rational interpolation, the performance of the approach presented here may be significantly degraded by the presence of noise in the available samples. This problem was not covered by this paper and is left for future investigations.