Structure-Preserving Interpolation of Bilinear Control Systems

In this paper, we extend the structure-preserving interpolatory model reduction framework, originally developed for linear systems, to structured bilinear control systems. Specifically, we give explicit construction formulae for the model reduction bases to satisfy different types of interpolation conditions. First, we establish the analysis for transfer function interpolation for single-input single-output structured bilinear systems. Then, we extend these results to the case of multi-input multi-output structured bilinear systems by matrix interpolation. The effectiveness of our structure-preserving approach is illustrated by means of various numerical examples.


Introduction
The modeling of various real-world applications and processes results in dynamical control systems usually including nonlinearities. Since linear approximations are very often incapable of capturing all the features of nonlinear systems, they are an insufficient description for use in optimization and controller design. A special class of nonlinear systems are bilinear control systems, which contain the multiplication of control and state variables, i.e., they are linear in state and control separately, but not together [28]. In the last decades, the class of bilinear systems became an essential tool in systems theory. They naturally appear in the modeling process of many physical phenomena, e.g., in the modeling of population, economical, thermal and mechanical dynamics [28,29], of electrical circuits [2], of plasma devices [30,31], or of medical processes [34]. Bilinear systems can also result from approximation of general nonlinear systems employing the Carleman linearization process [16,26]. Moreover, bilinear systems are nowadays often used in the parameter control of partial differential equations (PDEs) [24,25]. Looking back to the linear case, bilinear systems can be used as a generalizing framework in the modeling of linear stochastic [11] and parametervarying systems [7,10,15], allowing the application of established system-theoretic tools such as model order reduction for those system classes.
In this paper, we focus on structured bilinear systems. Those structures arise from the underlying physical phenomena. For example, in case of bilinear mechanical systems, one has the bilinear control system defined by Mq(t) + Dq(t) + Kq(t) = m j=1 N p,j q(t)u j (t) + m j=1 N v,jq (t)u j (t) + B u u(t), with M, D, K, N p,j , N v,j ∈ R n×n for all j = 1, . . . , m, B u ∈ R n×m and C p , C v ∈ R p×n . In (1), q(t) ∈ R n , u(t) ∈ R m , and y(t) ∈ R p denote, respectively, the states (degrees of freedom), inputs (forcing terms), and the outputs (quantities of interest) of the underlying dynamical system. Due to the usual demand for increasing accuracy in applications, the number of differential equations n, describing the dynamics of systems as in (1) quickly increases, resulting in a high demand on computational resources such as time and memory. One remedy is model order reduction: a new, reduced, system is created, consisting of a significantly smaller number of differential equations than needed to define the original one while still accurately approximating the input-to-output behavior. Then one can use this lower-order approximation as a surrogate model for faster simulations or the design of controllers. The classical (unstructured) bilinear first-order systems are described by the state-space representation with E, A, N j ∈ R n×n for all j = 1, . . . , m, B ∈ R n×m and C ∈ R p×n . There are different methodologies for model reduction of (2), e.g., the bilinear balanced truncation method [2,11,23], different types of moment matching approaches for the underlying multi-variate transfer functions in the frequency domain [3,5,14,17,18], the interpolation of complete Volterra series [8,19,35] or even the construction of reduced-order bilinear systems from frequency data with the bilinear Loewner framework [4,21].
While it is possible to rewrite (1) into a classical bilinear system (2), the original structure is completely lost, which can lead to undesirable results in terms of accuracy, stability or physical interpretation. Moreover, some other structured bilinear systems, such as those with internal delays (see Section 2.2.2), cannot be represented in the form (2). Therefore, here we develop a structure-preserving model reduction approach for different system structures involving bilinear terms. Following [6], which studied structured linear dynamical systems, our goal is to generalize the structured interpolation approach to a general set of multivariate transfer functions associated with different structured bilinear control systems to preserve the system structure in the reduced-order model. The question we aim to answer is how we can construct an interpolatory reduced-order model of, e.g., (1), that has the same structure. Towards this goal, we develop a structure-preserving interpolation framework for this special class of nonlinear systems, namely the structured bilinear control systems; thus extending the theoretical analysis and computational framework developed by [6] for linear systems to bilinear control systems.
In Section 2, we review the theory for classical first-order bilinear systems and motivate the more general structure, we will consider, via two examples. Section 3 gives subspace construction formulae for interpolatory model reduction bases in the case of single-input single-output (SISO) systems and illustrates the effectiveness of the approach employing two numerical examples. The developed theory is then extended further in Section 4 to the multi-input multi-output (MIMO) case by matrix interpolation. Section 5 concludes the paper.

Structured bilinear systems
In this section, we present the basic properties of the structured bilinear systems considered in this paper. To clarify the presentation, we first revisit the unstructured (classical) bilinear control systems as given in (2) and then generalize these concepts to the structured case.

Revisiting the classical first-order bilinear systems
Given the unstructured bilinear system (2), define N = N 1 . . . N m and let I m k be the identity matrix of dimension m k . Assuming for simplicity E to be invertible, the initial condition x(0) = 0, and some additional mild conditions, the output of (2) can be expressed in terms of a Volterra series [33], i.e., where g k , for k ≥ 1, is the k-th regular Volterra kernel given by Using the multivariate Laplace transform [33], the regular Volterra kernels (3) yield a representation of (2) in the frequency domain by the so-called multivariate regular transfer functions with s 1 , . . . , s k ∈ C. This compact expression is actually the collection of the different combinations of the bilinear matrices, i.e., we can write (4) as For SISO systems, (4) simplifies to As stated in Section 1, for the unstructured bilinear system case (2), there are already different model reduction techniques. For the structured bilinear systems we consider here, we will concentrate on interpolatory methods.
Note that the assumption of E being invertible is only made for ease of presentation. The interpolation theory and interpolatory properties of the reduced-order model developed in the following sections hold for the general situation, yet the final construction of the reduced-order model might need some additional treatment as in the linear and unstructured bilinear cases; see, e.g., [1,12,22].

Moving from classical to structured bilinear systems
For the transition from unstructured to structured bilinear systems, we start by recalling the case of linear systems. The classical (unstructured) linear dynamical systems are described, in state-space, by with E, A ∈ R n×n , B ∈ R n×m and C ∈ R p×n . Assuming the initial condition Ex(0) = 0, the Laplace transform maps this problem to the frequency domain: where X(s), U (s) and Y (s) denote the Laplace transforms of the time-dependent functions x(t), u(t), and y(t). Inspired by much richer structured systems than (6) appearing in the linear case such as those describing the dynamic response of a viscoelastic body, [6] introduced a more general system of equations in the frequency domain, given by with matrix-valued functions K : C → C n×n , B : C → C n×m and C : C → C p×n . Note that (7) contains (6) as a special case. Assuming the problem to be regular, i.e., there exists an s ∈ C for which the matrix functions are defined and K(s) is full-rank, the problem (7) leads to the general formulation of structured transfer functions of linear systems describing the input-to-output behavior in the frequency domain. Inspired by (8) for k ≥ 1 and where N (s) = N 1 (s) . . . N m (s) with the matrix functions C : C → C p×n , K : C → C n×n , B : C → C n×m , N j : C → C n×n for j = 1, . . . , m. This general formulation includes transfer functions of classical bilinear systems (4) since we can choose For the construction of structured reduced-order bilinear models, we will use the projection approach, i.e., we will construct two model reduction bases W, V ∈ C n×r such that the reduced-order bilinear system quantities will be computed by for j = 1, . . . , m. The structured reduced-order bilinear control system G is then given by the underlying reduced-order matrices from (10) and with the corresponding structured multivariate subsystem transfer functions for k ≥ 1.

Bilinear second-order systems
We revisit the example of second-order bilinear systems (1) given in Section 1. First, we note that (1) can be rewritten in the first-order (unstructured) form (2) by introducing the new state vector x(t) = [q T (t),q T ] T such that we obtain for any invertible matrix J ∈ R n×n . For this first-order companion realization (11), we know the frequency domain representation to be given by the multivariate regular transfer functions (4). If we now plug in the structured matrices from (11), we can make use of those special block structures.
In general, we obtain for the frequency-dependent center terms and, therefore, Using this, we obtain for the first part of the k-th regular transfer function where we used the notion N p = N p,1 . . . N p,m and N v = N v,1 . . . N v,m . Multiplication with the remaining terms yields the regular transfer functions of (1) to be written in the form Having the general formulation of regular transfer functions (9) in mind, we see that we can rewrite (12) in the structured bilinear form (9) by setting Now assume that we construct model reduction bases W and V and compute the reduced order model by projection as in (10). This leads to the reduced-order bilinear system Note that the reduced-order bilinear system in (13) has the same structure as the original one and can be viewed as a reduced second-order bilinear system, where the full-order matrices in (1) are simply replaced by the reduced analogues from (13).

Bilinear time-delay systems
Another structured bilinear control system example is the case of bilinear systems with an internal time-delay, i.e., for a delay 0 ≤ τ ∈ R, which has regular transfer functions of the form see [21]. As in the case of the bilinear second-order systems, we see that (14) can be written in the setting of (9) using As in Section 2.2.1, once the model reduction bases W and V are constructed, the resulting reduced-order model retains the delay structure of the original system as it is given by In Sections 3 and 4, we will show how to construct the model reduction bases W and V such that the structured reduced-order bilinear control system provides interpolation of the full-order subsystems.

Interpolation of single-input single-output systems
In this section, we assume the SISO system case, i.e., m = p = 1. Therefore, the bilinear part consists of, at most, one term N = N 1 and the matrix functionals C and B map frequency points only onto row and column vectors, respectively. In this setting, the regular transfer functions drastically simplify since (9) can now be written as for k ≥ 1. In the remainder of this section, we develop the theory for structure-preserving interpolation (both the case of simple and high-order (Hermite) interpolation) and then present numerical examples to illustrate the analysis.

Structured transfer function interpolation
We want to construct the model reduction bases W and V and the corresponding reduced structured-bilinear system via projection as in (10) such that its leading regular transfer functions interpolate those of the original one; i.e., G k (σ 1 , . . . , σ k ) = G k (σ 1 , . . . , σ k ), where σ 1 , . . . , σ k ∈ C are some selected interpolation points.
The following two theorems answer the question of how the model reduction bases V and W can be constructed independent of each other. In other words, the interpolation conditions are satisfied only via V or W , no matter how the respective other matrix is chosen. First, we consider the model reduction basis V .
Theorem 1 (Interpolation via V ). Let G be a bilinear SISO system, described by (15), and G the reduced-order bilinear SISO system constructed by (10). Let σ 1 , . . . , σ k ∈ C be interpolation points for which the matrix functions C(s), K(s) −1 , N (s) and B(s) are defined and K(s) is full-rank.
and let W be an arbitrary full-rank truncation matrix of appropriate dimension. Then the subsystem transfer functions of G interpolate those of G in the following way: Proof. First, we note that the constructed vectors are given by and that by construction all those vectors are contained in span(V ). Therefore, for the first transfer function we obtain where we used the fact that P v1 is an oblique projector onto span(V ), i.e., z = P v1 z holds for all z ∈ span(V ), and that K( Considering the second transfer function, we get using the same arguments as for the first transfer function and additionally the construction of v 2 and the oblique projector P v2 . Continuing with this argumentation, the desired result follows by induction over the transfer function index k. The proof of Theorem 1 shows that the recursive construction of the truncation matrix is necessary for the interpolation of higher-order transfer functions. Also, it should be noted that W was an arbitrary full-rank truncation matrix of suitable dimensions but with no additional constraints for the interpolation of (15). Theorem 2 is the counterpart to Theorem 1 by only giving constraints for the left model reduction basis W , while V is now allowed to be arbitrary.
Theorem 2 (Interpolation via W ). Let G, G, and the interpolation points σ 1 , . . . , σ k ∈ C be as in Theorem 1. Construct W using and let V be an arbitrary full-rank truncation matrix of appropriate dimension. Then the transfer functions of G interpolate the transfer functions of G in the following way: Proof. The proof of this theorem follows analogous to the proof of Theorem 1. We only need to note that the left projection space span(W ) involves the C(s) matrix, which takes always the last argument of G k into account. Therefore, the order of the interpolation points is reversed and the recursion formula follows the transfer function order going from left to right. The rest follows as in the proof of Theorem 1 by taking the Hermitian conjugate of the matrix functions for the construction.
The main difference between Theorem 1 and Theorem 2 is the order in which the interpolation points have to be used. Switching between the two projection schemes leads to a reverse ordering of the interpolation points for the intermediate transfer functions.
The last theorem of this section states now the combination of Theorem 1 and Theorem 2 by two-sided projection.
Theorem 3 (Interpolation by two-sided projection). Let G and G be as in Theorem 1 and let V be constructed as in Theorem 1 for a given set of interpolation points σ 1 , . . . , σ k ∈ C and W as in Theorem 2 for another set of interpolation points ς 1 , . . . , ς θ ∈ C, for which the matrix functions C(s), K(s) −1 , N (s) and B(s) are defined and K(s) is full-rank. Then the transfer functions of G interpolate the transfer functions of G in the following way: and additionally, Proof. Since the interpolation conditions in (16) follow directly from Theorem 1 and Theorem 2, we only need to prove (17), the mixed interpolation conditions. For q and η as described in the theorem, we obtain where we used the construction of span(V ) in the third and of span(W ) in the fourth lines as denoted and following the strategy in the proof of Theorem 1.
It is an important observation that we can interpolate higher-order transfer functions by only evaluating lower ones for the construction of the model reduction bases. Following Theorem 3, we can in fact interpolate transfer functions up to order k + θ. Also, we recognize that the two-sided projection-based interpolation is able to match k + θ + k · θ interpolation conditions at the same time. Those results are similar to the unstructured systems case [3]. The special case of identical sets of interpolation points is discussed in the following section regarding Hermite interpolation.

Hermite interpolation
As in the linear case, we can use the projection framework to interpolate not only the transfer functions but also their derivatives. In the setting of the multivariate transfer function appearing in bilinear systems, this amounts to partial derivatives with respect to the different frequency arguments. For ease of notation, we introduce an abbreviation for partial derivatives denoting the differentiation of an analytic function f : C k → C with respect to the variables s 1 , . . . , s k and evaluated at z 1 , . . . , z k ∈ C. Moreover, the Jacobian of f is denoted by as the concatenation of all partial derivatives.
The following theorem states a Hermite interpolation result via V only.
Theorem 4 (Hermite interpolation via V ). Let G be a bilinear SISO system, described by (15), and G the reduced-order bilinear SISO system constructed by (10). Let σ 1 , . . . , σ k ∈ C be the interpolation points for which the matrix functions C(s), K(s) −1 , N (s) and B(s) are analytic and and let W be an arbitrary full-rank truncation matrix of appropriate dimension. Then the transfer functions of G interpolate the transfer functions of G in the following way: Proof. First, we note that the case k = 1 was already proven in [6] and 1 = . . . = k = 0 corresponds to Theorem 1. For k = 2, we start with j 2 = 0 to investigate the partial derivative with respect to s 1 involving the bilinear term. Using the product rule, the partial derivative can be written as for some appropriate constants c i1 , c i2 ∈ C, and i 1 , i 2 = 0, . . . , 1 . Now, we can show where we first used the construction of v 1,j1 and then that of v 2,0 with the projector P v2,0 onto span(V ). By induction over j 2 , the results for the case k = 2 follow from [6]; and by induction over k and j k , using the same arguments, the rest of the theorem follows.
We note the difference between Theorem 1 and Theorem 4 in terms of the subspace construction. While for the previous interpolation results, we are able to recursively construct the next part of the model reduction subspace by using the previous one, this is not possible in Theorem 4 due to the frequency dependence of the bilinear term N (s). Also, it follows that for the interpolation of the -th derivative, = 1 + . . . + k , of the k-th transfer function G k in the interpolation points σ 1 , . . . , σ k , the minimal dimension of the projection space span(V ) is given by + k.
As before, we can consider the counterpart to Theorem 4. In addition to reversing the order of interpolation points, the order of the derivatives needs to be reverted as well for the Hermite interpolation.
Theorem 5 (Hermite interpolation via W ). Let G, G the original and reduced-order models, respectively, and the interpolation points σ 1 , . . . , σ k ∈ C be as in Theorem 4. Construct W using and let V be an arbitrary full-rank truncation matrix of appropriate dimension. Then the transfer functions of G interpolate the transfer functions of G in the following way . . .
Proof. Observing that the order of the derivatives changed in the same way as the interpolation points, the proof works analogously to the proof of Theorem 4 while building on the ideas from the proof of Theorem 2.
An interesting fact in the structured linear case, as stated in [6], is the implicit matching of Hermite interpolation conditions without sampling the derivatives of the transfer function. Next, we extend this construction to the structured bilinear case. This result becomes a special case of Theorem 3 by using identical sets of interpolation points for V and W .
Theorem 6 (Implicit Hermite interpolation by two-sided projection). Let G and G be as in Theorem 4. Also let V and W be constructed as in Theorems 1 and 2, respectively, for the same set of interpolation points σ 1 , . . . , σ k ∈ C, for which the matrix functions C(s), K(s) −1 , N (s) and B(s) are analytic and K(s) is full-rank. Then the transfer functions of G interpolate the transfer functions of G in the following way: and additionally, Proof. While most of the results directly follow from Theorem 3 by using identical sets of interpolation points for V and W , the Hermite interpolation of the complete Jacobian ∇G k of the k-th order transfer function is new. Since k = 1 (the linear subsystem) is covered by [6], we assume k > 1. Therefore, and by the structure of the multivariate transfer functions G k , three different cases can occur depending on the differentiation variable, i.e., we have as possible derivative terms. Since those three cases work analogously to each other, we restrict ourselves, for the sake of compactness, to the first one. First, we extend the expression of the partial derivative further into Therefore, for the complete partial derivative, we obtain where we used, as denoted by the underbraces, the construction of either span(W ) or span(V ), and the fact that the model reduction bases V and W are constant matrices. As stated before, the results for the other partial derivatives follow analogously, which proves interpolation of the full Jacobian in the end.
As in the previous section, by using two-sided projection we can match interpolation conditions for a larger number of interpolation points and higher-order transfer functions. Following the results of Theorem 3 we can expect, using derivatives for the two-sided projection, to match at least (k + ) + (θ + ν) + (k + ) · (θ + ν) transfer function values, where k, relate to span(V ) and θ, ν to span(W ), and where = 1 + . . . + k and ν = ν 1 + . . . + ν θ denote the orders of the partial derivatives and k, θ the orders of the transfer functions to interpolate. two-sided projection). Let G and G be as in Theorem 4 and let V be constructed as in Theorem 4 for a given set of interpolation points σ 1 , . . . , σ k ∈ C and orders of partial derivatives 1 , . . . , k , and W as in Theorem 5 for another set of interpolation points ς 1 , . . . , ς θ ∈ C and orders of partial derivatives ν 1 , . . . , ν θ , for which the matrix functions C(s), K(s) −1 , N (s) and B(s) are analytic and K(s) has full-rank. Then the transfer functions of G interpolate the transfer functions of G in the following way:
For an easier understanding of Theorem 7, we consider here a small theoretical example, where we only interpolate the linear part choosing k = θ = 1, the interpolation points σ, ς and for the partial derivatives = 1 = 2 and ν = ν 1 = 1. Then using the first part of Theorem 7 we enforce interpolation of the following terms by means of span(V ): , And similarly via span(W ), we enforce interpolation of By using two-sided projection, we can now additionally match higher-order transfer functions and their partial derivatives, namely G 2 (σ, ς), ∂ s1 G 2 (σ, ς), ∂ s2 G 2 (σ, ς), ∂ s 2 1 G 2 (σ, ς), ∂ s1s2 G 2 (σ, ς), ∂ s 2 1 s2 G 2 (σ, ς). As already realized in Theorem 6, two-sided projection with the same sets of interpolation points leads to additional interpolation of derivatives. This also works in combination with Theorem 7. The following corollary states a particular special case. Corollary 1. Assume G and G are constructed as in Theorem 7 for identical sets of interpolation points σ 1 , . . . , σ k ∈ C and matching orders of the partial derivatives, i.e., 1 = ν 1 , . . . , k = ν k . Then additionally to the interpolation results of Theorem 7 it holds Proof. The proof follows directly from Theorem 6 by setting the last partial derivative as the final interpolation condition of the left and right projection spaces.

Numerical examples
We illustrate the SISO analysis using two numerical examples, having the structured bilinearities as in Sections 2.2.1 and 2.2.2. We compare our resulting structure-preserving interpolation framework to other approaches from the literature that have been used to approximate structured bilinear systems without preserving the structure, as in, e.g., [2,21].
We compare the approximation error both in time and frequency domains. In time domain, we display a point-wise relative output error for a given input signal, namely for t ∈ [0, t f ], and in frequency domain, we display the point-wise relative error of the first and second subsystem transfer functions, i.e., The experiments reported here have been executed on a machine with 2 Intel(R) Xeon(R) Silver 4110 CPU processors running at 2.10 GHz and equipped with 192 GB total main memory. The computer runs on CentOS Linux release 7.5.1804 (Core) using MATLAB 9.7.0.1190202 (R2019b).

Damped mass-spring system
First, we consider a damped mass-spring system. The linear parts of the dynamics are modeled as in [27], describing a chain of masses connected by springs and dampers, where each mass is additionally connected to a separate spring and damper. In order to focus on only the mechanical structure, we removed the holonomic constraint from [27]. For the bilinear part, the springs are modeled to be dependent on the applied external force, such that a displacement to the right increases the stiffness due to compression of the springs and to the left decreases it due to the appearing strain. This results in a structured bilinear control system of the form with M, D, K, N p ∈ R n×n and B u , C T p ∈ R n . The input matrix is chosen to apply the external force only to the first mass, i.e., B = e 1 , and the output gives the displacement of the second mass, i.e., C = e T 2 , where e i denotes the i-th column of the identity matrix I n . The bilinear term is a scaled version of the stiffness matrix where S is a diagonal matrix containing entries as linspace(0.2, 0, n). For our experiment, we have chosen the original system to consist of n = 1 000 masses.
We construct three reduced-order models: (i) our structure-preserving bilinear interpolation, denoted by StrInt, (ii) two unstructured classical bilinear approximations by converting (18) to first-order form (11) followed by interpolatory model reduction of this first-order system, denoted by FOInt. Note that FOInt yields a reduced-order model of the form (11)  the underlying physical structure. Also, it needs to be remarked that the computational effort for the construction of FOInt is higher than for the structure-preserving approach due to solving underlying linear systems of doubled size, even in a structure exploiting implementation; see, e.g., [13]. Since the original system is a mechanical model, we use only a one-sided projection to preserve the mechanical properties in the reduced-order model, i.e., we apply Theorem 1 and set W = V . For all approximants, we focus on the first and second transfer functions and choose purely imaginary interpolation points. We construct StrInt and FOInt(12) by using the interpolation points ±logspace(-4, 4, 3)i such that the resulting reduced-order bilinear systems are of order r = 12, giving two different interpolations in the same frequency points. Since bilinear secondorder systems can be rewritten as first-order systems by doubling the state-space dimension, we construct additionally a second unstructured approximation FOInt(24) of order r = 24 by using ±logspace(-4, 4, 6)i, which has twice the order of StrInt. Figure 1 shows the time output of the original system, as well as that of the structure-preserving (StrInt) and first-order interpolations (FOInt(12), FOInt(24)), where we applied the input signal u(t) = sin(200t) + 200, which can be seen as a step signal with a sinusoidal disturbance. We see that while all three outputs are indistinguishable in the beginning, FOInt(12) becomes unstable after approximately 20 time steps and FOInt(24) after around 50 time steps, while StrInt accurately approximates the original system over the whole time range of interest. Even though the linear dynamics in FOInt (12) and FOInt (24) are asymptotically stable, these reduced-order models completely lack the underlying physical mechanical structure and they become unstable for the chosen input signal. On the other hand, by using one-sided projection, StrInt preserves all the mechanical (and physical) properties    of the original system in terms of symmetry and definiteness of the system matrices, which then leads to the stable time simulation behavior in this case. Figures 2 and 3 show the approximation results in the frequency domain for the first two transfer functions. Comparing StrInt and FOInt (12), the structure-preserving approximation is orders of magnitude better than the unstructured approximation of the same size. StrInt and FOInt(24) behave mainly the same, while, for higher frequencies, we can observe a numerical drift-off of the unstructured approximation.

Time-delayed heated rod
This example, taken from [21], models a semi-discretized heated rod with distributed control and homogeneous Dirichlet boundary conditions, which is cooled by a delayed feedback and is described by the PDE with (ζ, t) ∈ (0, π) × (0, t f ) and boundary conditions v(0, t) = v(π, t) = 0 for t ∈ [0, t f ]. After a spatial discretization using central finite differences, we obtain a bilinear time-delay system of the formẋ with A, A d , N ∈ R n×n , B, C T ∈ R n , and where we have chosen n = 5 000 for our experiments.  To compare with our structure preserving approximation (StrInt), in this example, we use the approach from [21] to construct an unstructured bilinear system (2) without time-delay using the bilinear Loewner framework, denoted by BiLoewner. For the structured interpolation, we have used the interpolation points ±logspace(-4, 4, 2)i for the first transfer function and ±logspace(-2, 2, 2)i for the second transfer function with the two-sided projection approach from Theorem 3. The resulting reduced-order bilinear time-delay system has order r = 8. For the bilinear Loewner method, we have chosen the interpolation points ±logspace(-4, 4, 80)i and used the rank truncation idea to obtain a classical (unstructured) bilinear system, also of order 8.
With the input signal u(t) = cos(10t) 20 + cos(5t) 20 , Figure 4 shows that (a) the output trajectories of the original system, the structure-preserving interpolation and the bilinear system without time-delay are indistinguishable in the eye ball norm (b) but the relative error reveals that StrInt is several orders of magnitude better than BiLoewner while having the same state-space dimension. The same behavior can be observed in the frequency domain for the first and second transfer functions as shown in Figures 5 and 6, i.e., by preserving the special structure of the original system we obtain a significantly better approximation of the same size.

Interpolation of multi-input multi-output systems
In this section, we will generalize the results from SISO structured bilinear systems to MIMO ones as in (9) and give a numerical example to illustrate the theory.

Matrix interpolation
In principle, all the results from Section 3 can directly be extended to the MIMO system case (9). However, one needs to realize that in this case, the quantities to be interpolated, i.e., the subsystem transfer functions, are matrix-valued. The main difference from the SISO case lies in the collection of the bilinear matrices into N (s) = N 1 (s) . . . N m (s) and the corresponding Kronecker products that produce the different combinations of the linear and bilinear parts in the k-th order transfer functions, e.g., in (5). Additionally, we will use the following notation as alternative way of concatenating the bilinear terms. In this paper, we will only focus on matrix interpolation, i.e., we will interpolate the full matrix-valued structured subsystem transfer functions. There is a concept of tangential interpolation [3,20] to handle matrix-valued functions in which interpolation is enforced only in selected directions. We will consider that framework in a separate work since the definition of tangential interpolation is not unified yet for bilinear systems [9,32], let alone the structured ones we consider here.
The following theorem extends the results from Theorems 1 to 3 to MIMO structured bilinear systems.
Theorem 8 (Matrix interpolation). Let G be a bilinear system, as described by (9), and G the reduced-order bilinear system, constructed by (10). Given sets of interpolation points σ 1 , . . . , σ k ∈ C and ς 1 , . . . , ς θ ∈ C, for which the matrix functions C(s), K(s) −1 , N (s), B(s) are defined and K(s) is full-rank, the following statements hold: (a) If V is constructed as then the following interpolation conditions hold true: . . .
(b) If W is constructed as then the following interpolation conditions hold true: (c) Let V be constructed as in part (a) and W as in (b), then, additionally to the results in (a) and (b), the interpolation conditions hold for 1 ≤ q ≤ k and 1 ≤ η ≤ θ.
Proof. Starting with part (a), we remember that the transfer functions can be rewritten by multiplying out the Kronecker products as For Hermite interpolation as in Theorems 4, 5 and 7, a similar extension to the MIMO case follows.
Proof. The results follow directly from Theorems 4, 5 and 7 with the same argumentation as in Theorem 8.
For completeness, also the implicit interpolation results are stated in the following corollary without additional proofs.
Corollary 2 (Two-sided matrix interpolation with identical point sets). Let G be a bilinear system, described by (9), and G the reduced-order bilinear system, constructed by (10). Given a set of interpolation points σ 1 , . . . , σ k ∈ C, for which the matrix functions C(s), K(s) −1 , N (s), B(s) are analytic and K(s) is full-rank, the following statements hold: (b) Let V and W be constructed as in Theorem 9 (a) and (b) for the interpolation points σ 1 , . . . , σ k , then additionally it holds

Numerical example
We illustrate the matrix interpolation results in a numerical example. The experiments reported here have been executed on the same machine and with the same MATLAB version as in Section 3.3. We reconsider the damped mass-spring system example from Section 3.3.1 with the following modifications: The mass, damping and stiffness matrices from (18) stay unchanged. The input forces are now applied to the first and last masses, i.e., the input term becomes B u = e 1 , −e n , and we observe the displacement of the second and fifth masses, which gives the output matrix C p = e 2 , e 5 T . Therefore, we have 2 inputs and outputs. We consider the same idea of bilinear springs as before but working in different directions, i.e., we have where S 1 is chosen, as before, as diagonal matrix with linspace(0.2, 0, n), and S 2 is chosen to be a diagonal matrix with linspace(0, 0.2, n) as entries. Overall, we have a damped massspring system of the form Mẍ(t) + Dẋ(t) + Kx(t) = N p,1 x(t)u 1 (t) + N p,2 x(t)u 2 (t) + B u u(t), with n = 1 000 masses for our experiments. As in Section 3.3.1, we compare the structure-preserving interpolation method (StrInt) with the unstructured one, using the first-order realization of (19) (FOInt). For the construction of StrInt and FOInt(36), we choose ±logspace(-4, 4, 3)i as interpolation points for the first transfer function and ±logspace(-3, 3, 3)i for the second one. Additionally, we construct another first-order approximation, FOInt(72), twice as large as the structured interpolation by taking ±logspace(-4, 4, 6)i and ±logspace(-3, 3, 6)i, as interpolation points for the first and second transfer functions, respectively. Also, we restrict ourselves again to a one-sided projection as in part (a) of Theorem 8 by setting W = V , which yields the reduced order r = 36 for StrInt and FOInt(36), and r = 72 for FOInt(72).
for t ∈ [0, 100]. The different lines in Figure 7a with the same color result from the two system outputs. In contrast to the SISO case, the linear part of the larger unstructured approximation (FOInt(72)) is not asymptotically stable anymore, which leads to the fast diverging behavior in the time simulation. The other first-order approximation (FOInt(36)) has a stable linear part but, as in the SISO case, is not able to produce stable results in the time simulation. StrInt again approximates the system's behavior accurately in the considered time range and, by using onesided projection, resembles the mechanical structures of the original system. Figures 8 and 9 show the results of the approximations for the first two transfer functions, where the relative errors are computed by for ω 1 , ω 2 ∈ [10 −2 , 10 2 ]. For both transfer function levels, we observe that FOInt(36) is not as accurate as StrInt and FOInt(72), which both nicely approximate the transfer functions except for higher frequencies, where the unstructured approximation seems to have the same numerical drift-off effect as in the SISO case.

Conclusions
We extended the structure-preserving interpolation framework to bilinear control systems. First, we developed the subspace conditions for structured interpolation for single-input single-output systems, both for simple and Hermite interpolation. These results were extended to structured multi-input multi-output bilinear systems as well in the setting of full matrix interpolation. The effectiveness of the proposed approach was illustrated for two structured bilinear dynamical systems: a mass-spring-damper system and a model with internal delay. The theory developed here can be applied to a much broader class of structures than these two examples.
In our examples, we made the rather simple choice of logarithmically equidistant interpolation points on the first two transfer function levels; thus the crucial problem of choosing good/optimal    interpolation points remains open. This question is not fully resolved even for structure-preserving interpolation of linear dynamical systems. Another issue to further investigate is the rapidlyenlarging reduced-order dimension in case of the matrix interpolation approach for multi-input multi-output systems. While in the linear case, tangential interpolation can be used to control the growth of the basis, there is no uniform treatment of tangential interpolation for bilinear systems yet. This issue will be studied in a separate work.