Structured interpolation for multivariate transfer functions of quadratic-bilinear systems

High-dimensional/high-fidelity nonlinear dynamical systems appear naturally when the goal is to accurately model real-world phenomena. Many physical properties are thereby encoded in the internal differential structure of these resulting large-scale nonlinear systems. The high-dimensionality of the dynamics causes computational bottlenecks, especially when these large-scale systems need to be simulated for a variety of situations such as different forcing terms. This motivates model reduction where the goal is to replace the full-order dynamics with accurate reduced-order surrogates. Interpolation-based model reduction has been proven to be an effective tool for the construction of cheap-to-evaluate surrogate models that preserve the internal structure in the case of weak nonlinearities. In this paper, we consider the construction of multivariate interpolants in frequency domain for structured quadratic-bilinear systems. We propose definitions for structured variants of the symmetric subsystem and generalized transfer functions of quadratic-bilinear systems and provide conditions for structure-preserving interpolation by projection. The theoretical results are illustrated using two numerical examples including the simulation of molecular dynamics in crystal structures.

Mathematics subject classification: 30E05, 34K17, 65D05, 93C10, 93A15 Novelty statement: We introduce new formulas for structured subsystem transfer functions to describe quadratic-bilinear systems with internal differential structures in the frequency domain.We formulate conditions on projection spaces to enforce structure-preserving interpolation for such transfer functions allowing for structure-preserving model reduction of quadraticbilinear systems.
Figure 1: Schematic illustration of the Toda lattice with particles [35].Atoms in a onedimensional crystal structure are represented as point masses and connected by exponential springs modeling the forces between the particles.

Introduction
The accurate modeling of real-world phenomena and processes yields dynamical systems typically including nonlinearities.Additionally, these systems often inherit some internal structure from the underlying physical nature of the considered problem.An example for such internal structures is the description of internal system states by second-order time derivatives as it is usually the case in the modeling of mechanical structures.Such nonlinear mechanical systems take the form M q(t) = f q(t), q(t), u(t) , with internal states q(t) ∈ R , describing the system behavior, the external controls u(t) ∈ R m that allow the user to change the internal behavior, and the quantities of interest y(t) ∈ R p that can be observed from the outside, e.g., by sensor measurements.Thereby, the first equation in ( 1) is a second-order differential equation with mass (descriptor) matrix M ∈ R × and the nonlinear time evolution function f : R × R × R m → R .The second, algebraic equation describes the quantities of interest as linear combination of the states and their first-order derivatives.Throughout this paper, we assume for any system to have homogeneous initial conditions.The Toda lattice is an example of a nonlinear mechanical system of the form (1). It is used in solid state mechanics to model the movement of particles in a one-dimensional crystal structure [33]; see Figure 1.The nonlinear time evolution function contains exponential terms that describe the forces between the different particles: f q(t), q(t), u(t) = − D q(t) −      e k 1 (q 1 (t)−q 2 (t)) − 1 e k 2 (q 2 (t)−q 3 (t)) − e k 1 (q 1 (t)−q 2 (t)) . . .e k q (t) − e k −1 (q −1 (t)−q (t)) with the positive semidefinite diagonal damping matrix D ∈ R × and the positive stiffness coefficients k 1 , . . ., k .See [35,Sec. 1.3.3]for the derivation of the differential model from the underlying Hamiltonian.In many applications, in particular those involving discretizations of partial differential equations, the number of differential equations in (1), describing the internal system behavior, is large and increases further with the demand for more accuracy.However, an increasing amount of differential equations also leads to an increasing demand for computational resources such as time and memory for simulations of the models or their use in optimization.A remedy to this problem is model order reduction, which aims for the computation of cheap-to-evaluate surrogate models described by significantly less differential equations, r , which approximate the input-to-output behavior of the original system as y − ŷ ≤ τ • u , in some suitable norms for the output of the reduced model ŷ and all admissible inputs u.
An established approach for model reduction of general (structured) nonlinear systems such as (1) is proper orthogonal decomposition (POD) [20,23,37], in which time simulations are used to extract information about the system dynamics for the construction of a basis matrix to project the system states.Other approaches aim for the extension of balancingrelated model reduction to nonlinear systems using Gramians defined via time simulations as in the empirical Gramian method [19,24,25] or by new energy measures [21,32] to construct suitable projection matrices.Nevertheless, a problem arising in the general system case is the approximation of the nonlinear time evolution function f in (1) that circumvents the computationally expensive lifting and truncation of the low-dimensional state in every time step.A solution to that are hyperreduction techniques such as the (discrete) empirical interpolation method ((D)EIM) [4,12,13,27], which computes a selection operator to restrict the evaluation of f to its rows most contributing.This introduces another layer of approximations and needs explicit access to the implementation of the original time evolution function f .
An alternative to hyperreduction that gained significant popularity in model reduction in the last decade is quadratic-bilinearization [18]; in optimization also known as Mc-Cormick relaxation [26]: For f smooth enough, general nonlinear systems can be rewritten into quadratic-bilinear form by introducing auxiliary variables and differential-algebraic equations, which then can be reduced directly using the classical projection-based model reduction approach.In the case of (1), the corresponding quadratic-bilinear system retains the internal mechanical structure as 0 = M q(t) + D q(t) + Kq(t) + H vv q(t) ⊗ q(t) + H vp q(t) ⊗ q(t) + H pv q(t) ⊗ q(t) + H pp q(t) ⊗ q(t) where and ⊗ denotes the Kronecker product.Due to the introduction of auxiliary variables, we have that n ≥ , which appears counter-intuitive to the actual task of reducing the number of internal system states in model reduction.However, the new nonlinearity structure of (3) allows to apply well-established model reduction techniques without the necessity of the hyperreduction step for the nonlinearity.The Toda lattice model (Figure 1) can also be rewritten into (3).The process of quadratic-bilinearization and the derivation of several structured quadratic-bilinear formulations for the Toda lattice model are shown in [35,Sec. 6.2].Furthermore, we refer the reader to the numerical experiments in Section 5.3 for more details.
So far, the literature on the model reduction of quadratic-bilinear systems mainly considered the case of unstructured, first-order systems of the form with E, A, N j ∈ R n×n , for j = 1, . . ., m, H ∈ R n×n 2 , B ∈ R n×m and C ∈ R p×n .Model reduction methods developed for (4) among others include the interpolation of multivariate subsystem transfer functions [1,2,6,18], Volterra series interpolation [2,9], balanced truncation [7], learning models from frequency domain data via the Loewner framework [15] or learning models from time domain data by operator inference [28,29].The reformulation of general nonlinear systems into quadratic-bilinear form (4) has also been proven to be an effective strategy for classical nonlinear model reduction methods such as POD [22].
In this paper, we extend the idea of quadratic-bilinear subsystem interpolation to systems with additional internal differential structures such as in the mechanical case (3).We propose extensions for the definitions of the first three symmetric subsystem transfer functions and the first three generalized transfer functions to the structured system case and then present subspace conditions for structure-preserving interpolation of these transfer functions.The effectiveness of resulting model reduction methods based on this interpolation theory is illustrated on two different structured examples including the Toda lattice model above.Parts of the theoretical results presented here were derived in the course of writing the dissertation of the corresponding author [35].
The rest of the paper is organized as follows: In Section 2, we recap the ideas of Volterra series expansions and unstructured quadratic-bilinear systems in the frequency domain.We present the definitions of structured transfer functions of quadratic-bilinear systems in Section 3 and the results on structure-preserving interpolation in Section 4. We employ the interpolation results for model reduction of two structured numerical examples in Section 5.The paper is concluded in Section 6.

Mathematical preliminaries
In this section, we recap the concept of Volterra series for quadratic-bilinear systems and two resulting transfer function formulations.

Volterra series expansions
The Volterra series expansion allows to describe the solution of nonlinear dynamical systems as a series of solutions of coupled linear systems [31].A common approach to derive the Volterra series expansion is by variational analysis [18].Let a scaled input signal αu(t), with α > 0, be given for the quadratic-bilinear system (4) and assume the system state to have an analytic representation of the form with a sequence of states x k (t).Inserting ( 5) into (4) and extracting the terms corresponding to the same power of the coefficient α yields the states x k (t) to be described by cascaded subsystems, which are linear in their respective unknown state x k (t).For example, the first three resulting linear subsystems for (5) are given by Applying the variation-of-constants formula to the subsystems in (6) allows the description of the input-to-output behavior of (4) via its Volterra series expansion: where g k (t 1 , . . ., t k ) are the Volterra kernels of the corresponding representation.The kernels used in (7) are of the symmetric type [18,31].Applying the multivariate Laplace transformation [31] to (7) results in an equivalent description of the quadratic-bilinear system (4) in the frequency domain by multivariate transfer functions.

Subsystem transfer functions of quadratic-bilinear systems
In this work, we restrict ourselves to the presentation of the transfer functions corresponding to the first three coupled linear subsystems for brevity and practical relevance.General formulas for arbitrarily high levels of multivariate transfer functions for (4) have been developed in [35, Sec.2.3.2].

Symmetric subsystem transfer functions
The symmetric subsystem transfer functions are based on the symmetric Volterra kernels from [31]; cf.(7).Historically, this is the first transfer function type that has been investigated for the model reduction of (4) in [18].Here, the term "symmetric" refers to the fact that the transfer function is invariant with respect to the order of its arguments.The first symmetric subsystem transfer function corresponds to the linear part of (4) as it can be seen in ( 6) such that with s 1 ∈ C. Thereby, the term g sym,1 (s 1 ) ∈ C n×m is used in the following for notational convenience and denotes the input-to-state transition of the first subsystem.The second symmetric subsystem transfer function depends on two complex frequency arguments s 1 , s 2 ∈ C and is given by where the function describing the input-to-state transition on the right-hand side of ( 9) is given by + N I m ⊗ g sym,1 (s 1 ) + g sym,1 (s 2 ) , with g sym,1 from (8), I m denoting the m-dimensional identity matrix and the column concatenation of the bilinear terms as Last, the third symmetric subsystem transfer function is defined similarly to (9) by with s 1 , s 2 , s 3 ∈ C and the input-to-state transition given via + N I m ⊗ g sym,2 (s 1 , s 2 ) + g sym,2 (s 1 , s 3 ) + g sym,2 (s 2 , s 3 ) , using g sym,1 from (8) and g sym,2 from (9).

Generalized transfer functions
In contrast to the symmetric case, the generalized transfer functions do not directly correspond to a Volterra kernel representation.They have been introduced in [15] for the extension of the data-driven Loewner framework to quadratic-bilinear systems and are inspired by the regular transfer functions of bilinear systems, which only consist of products of the terms of the dynamical system; see, e.g., [3].The formulation given in [15] for the single-input/single-output (SISO) system case has been extended to multi-input/multioutput systems in [35,Sec. 2.3.2].
As in the symmetric case, the first regular transfer function corresponds to the linear system components, with Also the second regular transfer function is uniquely defined resulting from one multiplication with the bilinear terms: For the third level however, two different generalized transfer functions exist.The first one is identical to the third regular bilinear transfer function with see also [10].The second one involves the quadratic term and is given by Note that the index levels of the generalized and symmetric transfer functions do not coincide since the second symmetric subsystem transfer function does contain the quadratic term in contrast to the second level generalized transfer function, which has only the bilinear terms.

Structured transfer functions of quadratic-bilinear systems
In this section, we extend the transfer function formulations from Section 2 to the setting of structured quadratic-bilinear systems starting with the introductory example of quadratic-bilinear mechanical systems.Based on this motivation, we introduce the formulas for structured symmetric and generalized transfer functions before we consider the case of quadratic-bilinear time-delay systems as another example for internal differential structures at the end of this section.

Transfer functions for mechanical systems
In general, any system in second-order form (3) can be rewritten in first-order form (4) using the concatenated first-order state x(t) = q(t) T q(t) T T .The first-order system matrices are then, for example, given by for j = 1, . . ., m, with the quadratic term Thereby, the matrix blocks in (17) are n × n matrix slices of the second-order quadratic terms in (3), e.g., for H pp modeling the multiplication of the state with itself we have with H pp,j ∈ R n×n for all j = 1, . . ., n.Now, we exploit the block structures of the matrices in ( 16) and ( 17) to derive the symmetric and generalized transfer functions of (3).In both transfer function cases, the first level transfer functions correspond to the linear system case and it can be observed that for (3) it holds that where g sym,1 (s 1 ) denotes here the input-to-state transition of the second-order state q.
For the next two levels, we concentrate first on the symmetric transfer function case.Inserting ( 16) and ( 17) into (9) yields the second symmetric subsystem transfer function of (3) to be where g sym,2 denotes the input-to-state transition of the second subsystem, g sym,1 is the input-to-state transition from the first subsystem in (18), and the bilinear terms are concatenated as Similarly, inserting the block matrices ( 16) and ( 17) into the unstructured third symmetric subsystem transfer function (11) leads to the structured third symmetric subsystem transfer function representation for (3).Due to the complexity of the resulting formula, we will not write it out here explicitly but outline some of its features.The third subsystem transfer function of (4) has a similar structure to ( 11) and ( 19) with linear combinations of the quadratic and bilinear terms multiplied with the previous input-to-state transition terms g sym,1 and g sym,2 from ( 18) and (19).Due to the occurrence of the sum of the frequency arguments in the first term of (19), this translates into the frequency dependence of the second-order quadratic and bilinear terms such that terms of the forms and appear.For more details, we refer the reader to the next section, which contains the definitions of the transfer function formulas for general structures.
For the generalized transfer functions, one can observe that the second level transfer function resembles the corresponding bilinear regular transfer function, for which the structured system case has been developed in [10].The resulting transfer function for (3) is thereby given as where the bilinear terms are concatenated as in (20).Similarly, the purely bilinear third level generalized transfer function of ( 3) is given by On the other hand, the third level generalized transfer function of (3) involving the quadratic term can be derived by inserting ( 16) into (15), which yields Overall, and similar to the linear and bilinear cases, we can observe the occurrence of the same terms describing the linear, bilinear or quadratic dynamics in the different transfer functions by means of the given system matrices and the complex variables s 1 , s 2 , s 3 .This motivates the definitions for the general structured framework in the upcoming section.

Structured transfer function formulas for quadratic-bilinear systems
Before deriving structured quadratic-bilinear transfer functions, we briefly recall the ideas from [5], which considered the structured transfer functions for linear dynamical systems in which frequency-dependent equations are used for describing the dynamics.First, consider linear first-order (unstructured) systems with the time domain representation By taking the Laplace transform of ( 25), the dynamical system in (25) can be equivalently described in the frequency domain as with s ∈ C, where U (s), X(s), Y (s) are the Laplace transforms of the inputs u(t), states x(t), and outputs y(t), respectively.Now consider a linear dynamical system with a second-order structure with the time domain representation As in the unstructured case, taking the Laplace transform of ( 27) yields the representation in the frequency domain Observe that in both cases of ( 26) and ( 28), the system states in the frequency domain can be described as solution of frequency-dependent linear systems of equations of the form where the matrix-valued functions K : C → C n×n , B : C → C n×m and C : C → C p×m describe the linear dynamics, and the input and output behavior of the system.In particular, we recover ( 26) by setting Similarly, we recover (28) by setting We refer the reader to [5] for further examples of structured dynamics that fit into the general framework of (29).Then, for every s ∈ C for which K(s) in ( 29) is invertible, the transfer function of the underlying linear dynamical system is given by For the extension to structured bilinear systems in [10], a new frequency-dependent function N : C → C n×nm was introduced, modeling the effect of the bilinear terms, where with N j : C → C n×n for all j = 1, . . ., m.In this manuscript, we further extend the structured transfer function framework to the quadratic-bilinear case, which appears ubiquitously in prominent applications as we briefly discussed in Section 1. First, we consider the symmetric transfer function case.Inspired by ( 8), (9), and ( 11), from the unstructured first-order case, and ( 18), (19), and ( 21), we introduce the following definition for the structured symmetric subsystem transfer functions.

Definition 1. Given matrix-valued functions of the form
for which there exists an s ∈ C at which they can be evaluated and K(s) is invertible.The first three structured symmetric subsystem transfer functions are defined as where the input-to-state transitions are recursively given by Similarly, we give a definition for the structured variant of the generalized transfer functions inspired by the first-order case ( 12)-( 15), and second-order case ( 18) and ( 22)-( 24) in the following.

Definition 2. Given matrix-valued functions of the form
for which there exists an s ∈ C at which they can be evaluated and K(s) is invertible.The first three levels of structured generalized transfer functions are defined as In [8], a simplified variant of the generalized transfer functions ( 12)-( 15) has been extended to systems with polynomial nonlinearities.Based on the structured definitions above, these simplified generalized transfer functions have then been extended to the structured case in [17].
Note that in both Definitions 1 and 2, the new matrix-valued function H : C×C → C n×n 2 results from the quadratic terms in the time domain.
Both examples of internal system structures considered so far can be represented in the new structured transfer function framework.Unstructured first-order systems of the form (4) are given by where the bilinear terms are concatenated as in (10).For the second-order system of the form (3), the symmetric and generalized transfer functions given in Section 3.1 can be recovered from Definitions 1 and 2 using with the bilinear terms concatenated as in (20).For the definition of higher level structured transfer functions for quadratic-bilinear systems see [35,Sec. 6.3].

Another example structure: Quadratic-bilinear time-delay systems
Before we consider interpolation of the structured transfer functions in Definitions 1 and 2, we present another example for internal system structures that are covered by the new structured transfer function framework.Quadratic-bilinear systems with constant time delays in the linear dynamic components can be written as with the matrices A k ∈ R n×n describing the effect of state delayed by τ k ∈ R ≥0 , for all k = 1, . . ., , and the remaining system matrices as defined in (4).Following the variational analyses from (6), we observe that the time-delay structure only affects the terms with the linear dynamics.Therefore, the structured transfer functions for (31) are given by using the matrix-valued functions in Definitions 1 and 2. This is in accordance to the results for bilinear time-delay systems obtained in [10,16].

Structured transfer function interpolation
In this section, we present results on the construction of structured interpolants for the symmetric or generalized transfer functions from the previous section.

Structure-preserving model reduction via projection
For the construction of interpolating reduced-order models, we will use the projection approach in this work.Thereby, two constant basis matrices V, W ∈ C n×r are constructed, which allow the computation of the reduced-order quantities via multiplication with the original system matrices.Given the full-order matrix-valued functions C : that describe a structured quadratic-bilinear system in the frequency domain, reduced-order model quantities are computed by where W H := W T denotes the conjugate transpose of the matrix W .The Kronecker product in the multiplication with the concatenation of the bilinear terms in (32) boils down to the multiplication of each single bilinear term with the two basis matrices as Moreover, the Kronecker product of the basis matrix V for the reduction of the quadratic term in (32) can be implemented efficiently without explicitly forming V ⊗ V , using techniques from tensor algebra; see, e.g., [6,7,35].Model reduction by projection preserves internal structures by construction.Any matrix-valued function can be decomposed into frequency-affine form, e.g., in the case of the term describing the linear dynamics, it can be written as with n K ∈ N, some frequency-dependent scalar functions h j : C → C and constant matrices K j ∈ C n×n , for all j = 1, . . ., n K .The reduced-order matrix-valued function is then given by with the reduced-order constant matrices K j ∈ C r×r , for all j = 1, . . ., n K .The frequencydependent scalar functions in (34) are the same as in (33), i.e., the internal structure is preserved and the reduced-order matrices replace their high-dimensional counterparts from the original system to describe the reduced-order model.
To illustrate the computation of reduced-order quadratic-bilinear systems via projection, we consider the two motivational differential structures from the previous section.In the case of second-order quadratic-bilinear systems (3), reduced-order systems are computed as with the reduced quadratic terms given by Evaluating the matrix products yields the reduced-order matrices that represent the reduced-order system in the same structure as the original system (3).Similarly, for the quadratic-bilinear time-delay system (31), reduced-order systems are computed via As in the second-order system case, evaluating the matrix products allows to replace the original, high-dimensional system matrices in (31) by the reduced ones to describe the reduced-order system using the same structure.
The essential question of projection-based model order reduction is the construction of the basis matrices V and W .In the following, conditions are derived to enforce interpolation of the original symmetric or generalized transfer functions by the corresponding transfer functions given via the reduced matrix-valued functions (32).

Interpolating symmetric transfer functions
In this section, we consider the interpolation of the structured symmetric subsystem transfer functions from Definition 1.The following proposition states some first results for the general interpolation of the first two symmetric subsystem transfer functions at different frequency points.
Proposition 1 ([35,Cor. 6.3]).Let G be a quadratic-bilinear system, described by its symmetric subsystem transfer functions G sym,k from Definition 1, and G the reduced-order quadratic-bilinear system constructed by (32), with its reduced-order symmetric subsystem transfer functions G sym,k .Also, let σ 1 , σ 2 ∈ C be interpolation points such that the matrix-valued functions C, K, B, N , H and K(.) −1 are defined in these points and their sum.Construct the basis matrix V by and let W be an arbitrary full-rank matrix of appropriate dimensions.Then, the symmetric subsystem transfer functions of G interpolate those of G in the following way: Note that using [35, Thm.6.2], the basis matrix V in Proposition 1 can be extended such that the third symmetric subsystem transfer function is interpolated as well.In practice however, this result is barely used due to the exponential increase of terms to evaluate for the construction of the projection space and the corresponding computational complexity.Therefore, it is omitted here.
The next proposition considers a similar interpolation result as in Proposition 1 by setting conditions on the second basis matrix as well.
Proposition 2 ([35, Lem.6.4]).Given the same assumptions as in Proposition 1, let the matrices V 1,1 and V 1,2 be as in Proposition 1. Construct the two basis matrices such that and let V and W be of the same dimension.Then, the symmetric subsystem transfer functions of G interpolate those of G in the following way: The result of Proposition 2 allows to enforce the same and more interpolation conditions than in Proposition 1 in an implicit way using the second basis matrix.This reduces the computational complexity of the construction of the basis matrices, since no nonlinear terms are involved, and allows to match more interpolation conditions with smaller reduced-order models.
The choice of interpolation points for the different subsystem levels is crucial for the quality of the computed reduced-order model.Good or even optimal choices of interpolation points are currently unknown.However, an advantageous choice to minimize the amount of basis contributions necessary for the interpolation of higher level subsystem transfer functions in the symmetric case is σ 1 = σ 2 = σ 3 = σ.The following theorem states the interpolation conditions for this particular selection of interpolation points and also gives conditions for the interpolation of the third subsystem transfer function.
Theorem 1.Let G be a quadratic-bilinear system, described by its symmetric subsystem transfer functions G sym,k as in Definition 1, and G the reduced-order quadratic-bilinear system constructed by (32), with its reduced-order symmetric subsystem transfer functions G sym,k .Also, let σ ∈ C be an interpolation point such that the matrix-valued functions C, K, B, N , H and K(.) −1 are defined at σ as well as at 2σ and 3σ.Construct the matrices and Then, the following statements hold true: (a) If the basis matrix V is such that and W is full-rank and of the same dimension as V , then the symmetric transfer functions of G interpolate those of G in the following way: (b) If the basis matrices V and W are such that span(V ) ⊇ span (V 1 ) and span(W ) ⊇ span (W 1 ) , and have the same dimension, then the symmetric transfer functions of G interpolate those of G in the following way: (c) If the basis matrices V and W are such that and both of appropriate dimensions, then the symmetric transfer functions of G interpolate those of G in the following way: Using the projector hold due to the construction of the space spanned by the columns of V via the columns of V 1 and V 2 .Using the projector H onto the space spanned by the columns of W , it also holds that such that the final result of the theorem holds via = G sym,3 (σ, σ, σ).

Interpolating generalized transfer functions
In this section, we investigate conditions on the projection spaces for the interpolation of the structured generalized transfer functions from Definition 2. The following proposition can be seen as an analog to Proposition 1 for the generalized case.
Proposition 3 ([35,Thm. 6.13]).Let G be a quadratic-bilinear system, described by its generalized transfer functions G (.) gen,k from Definition 2, and G the reduced-order quadraticbilinear system constructed by (32), with its reduced-order generalized transfer functions G (.) gen,k .Also, let σ 1 , σ 2 , σ 3 ∈ C be interpolation points such that the matrix-valued functions C, K, B, N , H and K(.) −1 are defined at these points.Compute and construct the basis matrix V such that Let W be an arbitrary full-rank matrix of appropriate dimensions.Then, the generalized transfer functions of G interpolate those of G in the following way: Remark 1.As for the symmetric subsystem transfer functions, it is possible to reduce the dimensions of the constructed projection space by choosing suitable interpolation points and, also as in the symmetric subsystem transfer function case, this can be achieved by choosing σ 1 = σ 2 = σ 3 .However, only marginal savings in terms of basis contributions and computational costs can be achieved by this since in Proposition 3, only the matrix V 1,2 can be omitted.
Note that in Proposition 3, the basis contribution from V 3,1 results in the interpolation of the third generalized transfer function with two bilinear terms G (N,(N,(B))) gen,3 . It may happen that no interpolation conditions are imposed for this transfer function such that the subspace dimensions can be reduced by omitting V 3,1 .As in the case of symmetric subsystem transfer functions, the second basis matrix W can be used to reduce the minimal dimension of the constructed subspaces and to simplify the construction of the subspaces.These results are given in the following corollary.
Corollary 1.Let G be a quadratic-bilinear system, described by its generalized transfer functions G (.) gen,k from Definition 2, and G the reduced-order quadratic-bilinear system constructed by (32), with its reduced-order generalized transfer functions G (.) gen,k .Also, let σ 1 , σ 2 ∈ C be interpolation points such that the matrix-valued functions C, K, B, N , H and K(.) −1 are defined at these points.Let the basis matrices V and W be constructed by and are of the same dimension.Then, the generalized transfer functions of G interpolate those of G in the following way: Proof.The result follows directly from [35,Lem. 6.15] by restriction to two interpolation points.
Similar to Proposition 2, the result in Corollary 1 states that the interpolation of higher level transfer functions is possible in an implicit way without evaluating any of the nonlinear terms.Corollary 1 shows the version of the implicit interpolation result with the smallest achievable minimal subspace dimensions for the interpolation of the third level generalized transfer function with quadratic term.The choice of identical interpolation points will not further reduce the dimensions of the projection spaces, but allows to replace the interpolation of the first level generalized transfer function in two points by matching the transfer function value and its derivative in one point; see [5] and [35,Lem. 6.15].

Numerical experiments
Now, we employ the interpolation results from above for constructing structured reducedorder quadratic-bilinear systems in two numerical examples.The experiments were run on compute nodes of the Greene high-performance computing cluster of the New York University using 16 processing cores of the Intel Xeon Platinum 8268 24C 205W CPU at 2.90 GHz and 32 GB main memory.We used MATLAB 9.9.0.1467703 (R2020b) running on Red Hat Enterprise Linux release 8.4 (Ootpa).The source code, data and results of the numerical experiments are open source/open access and available at [36].

Experimental setup
In both numerical examples, we compute reduced-order models via structure-preserving interpolation of the symmetric subsystem transfer functions and the generalized transfer functions, denoted by SymInt and GenInt, respectively.We compute models using either (i) only the construction of the basis matrix V and a one-sided projection by setting W = V , which we abbreviate further on by V, or (ii) by also constructing the left basis matrix W for a two-sided projection following the results in Theorem 1 and Corollary 1, which we abbreviate by VW.For simplicity, the interpolation points are chosen logarithmically equidistant on the imaginary axis in all cases.If we compute only the basis matrices for interpolation without additional information, this is denoted by equi.On the other hand, if we oversampled the frequency range of interest and compressed the resulting basis to a prescribed dimension, e.g., using pivoted QR, this is denoted by avg; cf.[35,Rem. 3.3].In all cases, we focus on the interpolation of either (i) the first two symmetric subsystem transfer functions or (ii) the first two levels of the generalized transfer functions and the third level generalized transfer function containing the quadratic term.The following overview summarizes the considered interpolation methods: SymInt(V, equi) is the interpolation of symmetric subsystem transfer functions via onesided projection by constructing the basis matrix V .
SymInt(VW, equi) is the interpolation of symmetric subsystem transfer functions via twosided projection.Additional interpolation points are selected for the construction of W to match the dimension of V .
SymInt(V, avg) is the approximation of an interpolation basis for symmetric subsystem transfer functions using only samples for the construction of V and one-sided projection.
SymInt(VW, avg) is the approximation of left and right interpolation bases for symmetric subsystem transfer functions using samples for the construction of V and W and twosided projection.Additional interpolation points are selected for the construction of W to match the computational work to the construction of V .
GenInt(V, equi) is the interpolation of generalized transfer functions via one-sided projection by constructing the basis matrix V .Samples from the second and third level transfer functions are taken alternating.
GenInt(VW, equi) is the interpolation of generalized transfer functions via two-sided projection.For the construction of V , samples from the second and third level transfer functions are taken alternating.Additional interpolation points are selected for the construction of W to match the dimension of V .
GenInt(V, avg) is the approximation of an interpolation basis for generalized transfer functions using only samples for the construction of V and one-sided projection.At all interpolation points, second and third level samples are taken.
SymInt(VW, avg) is the approximation of left and right interpolation bases for generalized transfer functions using samples for the construction of V and W and two-sided projection.For the construction of V , second and third level samples are taken at all interpolation points.Additional interpolation points are selected for the construction of W to equalize the computational work to the construction of V .
As an additional comparison, we have computed reduced-order models via proper orthogonal decomposition (POD).All POD models have been trained via simulations of the unit step response to remove the correlation of the training and test input signals.For a fair comparison, the trajectory lengths used for POD are chosen with respect to comparable amounts of computational work to the interpolation methods.We consider therefore: POD, which has a computational workload similar to SymInt(V, equi) and GenInt(V, equi), and POD(avg), which uses similar to SymInt(V, avg) and GenInt(V, avg) an oversampling and computes the orthogonal basis via truncated singular value decomposition.
For the comparison of the reduced-order models in time domain, we simulate the models over finite time intervals using input signals taken from a Gaussian process GP(µ, K), with constant mean µ ∈ R and the squared exponential kernel where ς ≥ 0 is a smoothing parameter.The parameters µ and ς are chosen independently for the two examples and are given below.For visualization, we compute and plot the maximum pointwise relative errors relerr(t) := max j y j (t) − ŷj (t) y j (t) .
Also, we compute discretized approximations of the relative L 2 and L ∞ errors via where y h , ŷh ∈ R p×n h are the discretized output signals of the original and reduced-order model, respectively, in the time interval [0, t f ], and vec(.) is the vectorization operator.
In frequency domain, we consider the pointwise relative spectral norm errors defined as over the limited frequency intervals ω, ω 1 , ω 2 ∈ [ω min , ω max ].Additionally, we compute approximations to the relative L ∞ -norm errors for the first and second symmetric subsystem transfer functions via relerr using 500 logarithmically equidistant sampling points in the frequency interval of interest [ω min , ω max ].
For further details on the experimental setup, we refer the reader to the accompanying code package [36].

Quadratic time-delayed reaction-diffusion model
As first example, we consider the time-delayed heated rod with bilinear feedback from [11,16], to which we append a quadratic reaction term to obtain with (t, ζ) ∈ (0, t f ) × (0, π), boundary conditions ν(t, 0) = ν(t, π) = 0 for all t ∈ [0, t f ], and the constant time delay τ = 1.After spatial discretization with central finite differences, we obtain a quadratic-bilinear time-delay system of the form For our experiments, we have chosen n = 2 000, m = 2 and p = 2, where u 1 controls the temperature of the first third of the rod and u 2 the rest, and y 1 observes the temperature of the first half of the rod and y 2 of the second half.The system has zero initial conditions x(t) = 0 for all t ≤ 0. The reduced-order models are computed as explained in Section 5.1 with a reduced order of r = 24.In Table 1, we can see that all reduced-order models perform well for this example.However, the interpolation-based models provide smaller errors than those generated by POD, where the better POD model computed by POD(avg) performs mildly worse than the worst interpolation-based reduced-order models.Comparing the different interpolation approaches we can observe that mostly, the interpolation of symmetric The best reduced-order models from each generating approach are shown.All reduced-order models can recover the system behavior for the given input signal, but the interpolation-based reduced-order models perform around four orders of magnitude better in terms of accuracy than the model generated by POD(avg).
transfer functions performs better in the time simulations and in frequency domain than the interpolation of the generalized transfer functions.For the reduced-order models that provide exact interpolation (equi), the sampling of higher-order terms in the generalized transfer function needed to be restricted to match the reduced basis dimension.This restriction is removed in avg such that more information about the bilinear and quadratic terms can be obtained in sampling the generalized transfer functions compared to the symmetric transfer function setting, but this does only lead to a better reduced-order model when measured with relerr L∞ .Figure 2 shows the time simulation of the full-order model and the best performing reduced-order models from each method over the time interval [0, 30] s.We restricted Figure 2a to only the first output signal for clarity, but the pointwise relative errors in Figure 2b are computed over both output entries.For the input signals, we have chosen the mean µ = 2 and the parameter ς = 0.25.The interpolation-based methods clearly outperform POD(avg) by approximately four orders of magnitude.There is no significant difference between the errors of SymInt(VW, avg) and GenInt(VW, avg).
Similar results can be observed in frequency domain.The first symmetric subsystem transfer functions are shown in Figure 3, while Figure 4 illustrates the pointwise relative errors of the second symmetric subsystem transfer functions.In both cases, we only show the best performing methods in the frequency interval ω, ω 1 , ω 2 ∈ [10 −3 , 10 3 ] rad/s.As in the time domain, the interpolation-based methods outperform POD(avg) by several orders of magnitude in terms of accuracy.Further plots of the other methods in time and frequency domain can be found in the accompanying code package [36].

Particle motion in one-dimensional crystal structures
As second example, we consider the motion of particles in a one-dimensional crystal structure described by the Toda lattice model from the introduction; see Figure 1 The best reduced-order models from each generating approach are shown.All reduced-order models can recover the system behavior for the given input signal, but the interpolation-based reduced-order models perform around six orders of magnitude better in terms of accuracy than the model generated by POD(avg).
By differentiating (35) twice, the Toda lattice model can be written as a system of quadratic-bilinear ordinary differential equations of the form 0 = M q(t) + D q(t) + Kq(t) + H vv q(t) ⊗ q(t) + H pv q(t) ⊗ q(t) with the dimensions as in (3) and m = 1 input and p = 1 output.The exact parameterization of the matrices in (36) can be found in [35,Sec. 6.5].For our experiments, we use the same setup as in [35, Sec.6.5] with = 2 000 particles such that (36) has the order n = 4 000.It has been observed in [35] that the internal block structures of the matrices in (36) resulting from the original and auxiliary variables should be preserved for stability of the reduced-order models.Therefore, we follow the suggestion in [35,Sec. 6.5] and use the split congruence transformation approach [14,30,34].That is, given a basis matrix and similarly for a left basis matrix W .The extended basis matrices are then used for model reduction by projection.We apply this approach in our experiments to modify all projection basis matrices computed as described in Section 5.1.By construction, it holds that span(V ) ⊆ span( V ).Therefore, if V was constructed to satisfy any subspace conditions in Section 4 for interpolation, the basis matrix V also satisfies these conditions such that interpolation properties are preserved.For the comparison in our experiments, we have chosen the reduced order of all computed models to be 2r = 120.The results are shown in Table 2.The best performing model in terms of time simulation error is SymInt(V, equi) followed by its counterpart for the generalized transfer functions GenInt(V, equi).None of the models resulting from a two-sided projection has a stable time-domain simulation, which appears to be a consequence of loosing additional mechanical properties by V = W . Also none of the POD generated models performs stable in the time-domain simulation.Here, the large frequency domain errors indicate that the approximation quality is not sufficient to approximate the system behavior well enough.In frequency domain, we observe similarly to the previous numerical example that SymInt(VW, avg) and GenInt(VW, avg) perform best in terms of approximating the first symmetric subsystem transfer function.
The time-domain simulations of the full-and the best performing reduced-order models from each method are shown in Figure 5 in the time interval [0, 100] s.For the input signal, we have chosen the mean µ = 0 and the smoothing parameter ς = 2. POD(avg) performs visibly stable only until around 60 s, while the other two models follow the system behavior over the complete time interval.The pointwise relative errors in Figure 5b reveal POD(avg) to be more accurate than SymInt(V, equi) and GenInt(V, equi) for the first half of the time interval before it assumes the same error level as the other two methods and finally becomes unstable.SymInt(V, equi) and GenInt(V, equi) have overall a similar error behavior, with the errors of GenInt(V, equi) being mildly larger.On the other hand, in frequency domain, we can observe a similar behavior compared to the previous numerical example.The results for the same reduced-order models that performed best in time domain can be seen in Figures 6 and 7, with the frequency interval of interest ω, ω 1 , ω 2 ∈ [10 −3 , 10 3 ] rad/s.The POD generated models perform worst, with the exception of GenInt(VW, equi).The models based on oversampling generalized transfer functions perform better than those based on oversampling the symmetric transfer functions due to the additional information obtained from the nonlinear terms.However, in this example we can observe that the oversampling procedure may produce larger errors than exact interpolation.In particular, models computed via avg provide worse approximation errors for larger frequencies, where the transfer functions are converging to zero, while the models with equi preserve the system behavior due to the exact interpolation.Further plots of the other methods in time and frequency domain can be found in the accompanying code package [36].

Conclusions
We have extended the structure-preserving interpolation framework to quadratic-bilinear systems.Based on two motivating structured examples, we have introduced the structured variants of the symmetric subsystem and generalized transfer functions of quadraticbilinear systems.For both transfer function types, we provided subspace conditions enabling the computation of interpolating structured reduced-order models by projection.The theoretical findings are then used to compute structured reduced-order models in two numerical examples.The theory presented here can be applied to a much broader class of structures than those used here for illustrations.

FOM
SymInt(V, equi) GenInt(V, equi) POD(avg) Figure 5: Time simulation of the Toda lattice example: The best reduced-order models from each generating approach are shown.Only the interpolation-based reduced-order models recover the system behavior over the full time interval, while POD(avg) becomes unstable after about 60 s.
The numerical results suggest that the interpolation of symmetric transfer functions provides more accurate reduced-order models than the generalized transfer functions for the same reduced order.This is most certainly a consequence of the restriction to only products of system terms in the generalized transfer functions.However, we have seen that when the basis matrices are first constructed by oversampling and then compressed, the generalized transfer function interpolation framework provides more accurate reducedorder models for the same computational costs as for the symmetric transfer functions due to additional information obtained from the nonlinear terms.The authors of [17] extend on this oversampling idea and the definition of structured generalized transfer functions for quadratic-bilinear systems from this paper to propose structure-preserving model reduction for systems with polynomial nonlinearities based on the interpolation of generalized transfer functions with at most one nonlinear component.While this gives a first efficient approach to simulation-free model reduction for polynomial systems, there are many open questions left.One related to our observations in this work is the question whether exact interpolation of transfer functions based on Volterra kernels may perform better than an oversampling procedure.This needs a thorough investigation of interpolation conditions for polynomial systems, which we will address in some future work.
Another transfer function type for quadratic-bilinear systems are regular subsystem transfer functions.These have been omitted in this paper since for the choice of identical interpolation points, the projection spaces of regular and symmetric subsystem transfer functions coincide.The formulas for structured variants of regular transfer functions and results on interpolation conditions will be presented in a separate work.
For simplicity of exposition, we have restricted the numerical experiments to only logarithmically equidistant interpolation points on the imaginary axis.While such a procedure is often sufficient in practice, the question of good or even optimal interpolation points remains open and crucial for the success of such model reduction methods.This is still an unresolved issue even in the case of structured linear systems and needs further investigation in the future.

Proof.
Part (a) is an extension of Proposition 1 to the third symmetric subsystem transfer function with the special selection of interpolation points and follows immediately from [35, Thm.6.2].Part (b) follows directly from Proposition 2 such that only Part (c) is left to be proven.The first three interpolation conditions follow from previous results, therefore we concentrate on the third symmetric subsystem transfer function.Inserting the selection of interpolation points into Definition 1 yields

Figure 2 :
Figure 2: Time simulation of the time-delay example: The best reduced-order modelsfrom each generating approach are shown.All reduced-order models can recover the system behavior for the given input signal, but the interpolation-based reduced-order models perform around four orders of magnitude better in terms of accuracy than the model generated by POD(avg).

Figure 3 :
Figure 3: Sigma plots showing G(iω)) 2 of the first symmetric subsystem transfer function of the time-delay example:The best reduced-order models from each generating approach are shown.All reduced-order models can recover the system behavior for the given input signal, but the interpolation-based reduced-order models perform around six orders of magnitude better in terms of accuracy than the model generated by POD(avg).

Figure 4 :
Figure 4: Second symmetric subsystem transfer function relative approximation errors relerr(ω 1 , ω 2 ) of the time-delay example: The best reduced-order models from each generating approach are shown.The errors of both interpolation-based reduced-order models are at least four orders of magnitude better than those of the POD(avg) model.
Maximum pointwise relative errors.

Table 1 :
Errors computed as shown in Section 5.1 for the time-delay example: POD computes the worst performing reduced-order model, while POD(avg) is only as good as the worst interpolation-based models.Here, SymInt(VW, avg) performs best w.r.t.three out of the four error measures, while according to relerr

Table 2 :
Error table of the Toda lattice example: The errors are computed as shown in Section 5.1.Only the interpolation methods with one-sided projections provide reduced-order models that are stable in the time simulation.The models with ∞ error are unstable.SymInt(V, equi) outperforms GenInt(V, equi) by around a factor of 2.