While the underlying dynamical system acts as a black box, model identification tools are important for the reliability of the discovered models (i.e., stability, prediction). At the same time, these discovered models might have large dimension and hence are not suitable for fast numerical simulation. The (LF) is a direct data-driven interpolatory method able to identify and reduce models derived directly from measurements. For measured data in the frequency domain, the (LF) is well established for linear and non-linear systems (e.g., bilinear or quadratic-bilinear systems) see [5,18]. In the case of time domain data, the (LF) was already applied for approximating linear models [20,26,17]. The mathematical description of dynamical systems uses the concept of DAEs1 which can simulate the input-u(t) to output-y(t) relation as in Fig. 1, where the differential operator denoted as f and the algebraic operator denoted as z.

Σ : x(t) = f(x(t), u(t)), y(t) = z(x(t), u(t)) input y(t) output Figure 1: Mathematical formalism for evolutionary phenomena.

In this section, we will briefly present some important material from system theory starting from the linear case.

Consider SISO linear, time-invariant systems with n internal variables (called "states" whenever the matrix E is non-singular).

E x(t) = Ax(t) + Bu(t),

where E, A ∈ R n×n , B ∈ R n×1 , C ∈ R 1×n . In the sequel we will restrict our attention to invertible matrix E and with zero D-term (D = 0) in the state-output equation2. The resulting system (1) is:

x(t) = Ax(t) + Bu(t), y(t) = Cx(t), t ≥ 0.

(

The explicit solution with the convolution integral3 notation and the time domain linear kernel h(t) as the impulse response of the system can be written as:

By assuming zero initial conditions and performing a Laplace transform, we obtain the transfer function description:

where Y (s), U(s) stand for the input and the output in the frequency domain.

A large class of nonlinear systems can be described by means of the Volterra-Wiener approach in [30].

Other relevant works on nonlinear systems and nonlinear modelling/identification include Schetzen (1980), Chen and Billings (1989), Boyd and Chua (1985) et. al.

1DAEs: Differential Algebraic Equations. 2The state-output equation often is represented as y(t) = Cx(t) + Du(t).

The aim in this study is to identify and reduce special types of non-linear systems (s.a., bilinear) from time domain measurements. By knowing only the input and the simulated or measured output in the time domain as in Fig. 2, we will identify the hidden model which explains this phenomenon. In such situations where only snapshots are available, beyond the linear fit which is very well established to a series of methods, a non-linear fit of a special type will be able to capture the hidden non-linearity.

Σ : unknown input y(t) output Figure 2: The input-output mapping from the data-driven perspective with the unknown system Σ. Although specific structures of the unknown system can be assumed/inspired by the physical problem. For instance, if the underlying physical phenomenon is fluid flow inside a control volume, quadratic models should be constructed e.g., [18].

The input-output relationship for a wide class of non-linear systems [30] can be approximated by a Volterra series up to a sufficiently high N as:

where

as the n th order Volterra kernel.

Definition 2.1. The n th order generalized frequency response function (GFRF) is defined as:

which is the multidimensional Fourier4 transform of

By applying the inverse Fourier transform of the n th order GFRF, Eq. 6 can be written as:

The n th Volterra operator is defined:

so that it holds y n = V n (u, u, ..., u).

Remark 2.1. Homogeneity of the Volterra operator The map u(t) → y n (t) is homogeneous of degree n, that is, αu → α n y n . Each Volterra kernel h n (t) determines a symmetric multi-linear operator. This will allow to order the non-linear terms in such a way that larger powers of the amplitude (α n ) will eliminate the corresponding dynamics term. That is precisely the sense of approximating weakly non-linear systems with Volterra series (i.e., with 0 < α < 1, it holds α n → 0 as n → ∞).

4As the frequency s = jω 1 lies on the imaginary axis, the Laplace transform simplifies to Fourier transform.

Consider the excitation of a system with an input consisting of two complex exponentials as in Eq. 9.

Such inputs are typically used in chemical engineering applications as [28].

By using the above input in Eq. 5, we can derive the first Volterra term with n = 1 as:

Similarly, for the second term we can derive:

Remark 2.2. (Conjugate symmetry):

The input amplitude is A, the angular frequency is ω, the imaginary unit is j, the first order FRF5 is H 1 ( jω), and H n ( jω, ..., jω), for n ≥ 2, are the higher-order FRFs or GFRFs. Then, the n th Volterra term can be written as:

n ( jω)e jω p, q t , ω p,q = (pq)ω.

where the following brief notations have been used:

;jω, ...,jω q-times ), ω p,q = (pq)ω, n C q = n! q!(nq)! .

The n th harmonic in the time domain is:

The n th harmonic in the frequency domain by applying Fourier transform in Eq. 14 is the following:

where δ(•) is the Dirac delta distribution.

5FRF: Frequency Response Function and GFRF: Generalized FRF. Figure 3: An instance of the single-sided power spectrum with a singleton input with ω = 1 is depicted. The underlying system is non-linear and as a result higher harmonics appeared with a DC (Direct Current -non-periodic) term as well.

We start with an account of the Loewner framework (LF) in the linear case [2,6,24]. The (LF) as an interpolation method seeks reduced models whose transfer function matches that of the original system at selected interpolation points. An important attribute is that it provides a trade-off between accuracy of fit and complexity of the model. It constructs models from given frequency data in a straightforward manner. In the case of SISO systems, we have the rational scalar interpolation problem to solve.

Consider a given set of data (abstract frequencies) as:

We partition the data in two disjoint sets:

where

The objective is to find H(s) ∈ C such that:

The left data is rearranged compactly as:

while the right data as:

Interpolation points are determined by the problem or are selected to achieve given model reduction goals. For ways of choosing the interpolation grids and of partitioning the data into the left and right sets, we refer the reader to the recent survey [24].

Given a row array of complex numbers (µ j , v j ), j = 1, . . . , n, and a column array, (λ i , w i ), i = 1, . . . , n, (with λ i and the µ j mutually distinct) the associated Loewner matrix L and the shifted Loewner matrix L s are defined as:

Definition 3.1. If g is rational, i.e., g(s) = p(s) q(s) , for appropriate polynomials p, q, the McMillan degree or the complexity of g is deg g = max{deg(p), deg(q)}. Now, if w i = g(λ i ), and v j = g(µ j ), are samples of a rational function g, the main property of Loewner matrices asserts the following. Theorem 3.1. [2] Let L be as above. If k, q ≥ deg g, then rank L = deg g. In other words the rank of L encodes the complexity of the underlying rational function g. Furthermore, the same result holds for matrix-valued functions g.

is a minimal realization of an interpolant for the data, i.e., H(s) = W(L s -sL) -1 V, interpolates the data. Otherwise, as shown in [2], the problem in Eq. 16 has a solution provided that

for all s ∈ {µ i } ∪ {λ j }. Consider then the short SVDs:

Remark 3.1. r can be chosen as the numerical rank (as opposed to the exact rank) of the Loewner pencil.

Theorem 3.2. The quadruple ( Ẽ, Ã, B, C) of size r × r, r × r, r × 1, 1 × r, given by: For more details on the construction/identification of linear systems with the (LF), we refer the reader to [6,24] where both the SISO and MIMO cases are addressed together with other more technical aspects (e.g., how to impose the construction of real-valued models, etc.).

In recent years, projection-based Krylov methods have extensively been applied for model reduction of bilinear systems. We mention the following contributions [5,1,7,9,10,11,14,16,29] and the references within.

Bilinear systems are described by the following set of matrices; E,A,N,B) and characterized by the following equations:

where

In what follows, we restrict our analysis to systems with non-singular E matrices. Moreover, this matrix can be considered to be the identity matrix.

The properties of the growing exponential approach can be adapted readily to the problem of finding transfer functions descriptions from constant-parameter (stationary) state equations. Let us consider the bilinear model in Eq. 19 with zero initial conditions. A single oscillatory input with amplitude A < 1 is considered as in Eq. 9.

where a = A/2 and a ∈ (0, ) with 0 < < 1/2 and for all t ≥ 0. The steady state solution for the differential equation in Eq. 19 can be written as:

The symbol G p,q n 6 is the n th input to state frequency response containing p-times the frequency ω and q-times the frequency -ω. By substituting in Eq. 19 and collecting the terms of the same exponential (as the e jω m t ), we can derive the input to state frequency responses G for every n. By denoting the resolvent Φ( jω) = ( jωE -A) -1 ∈ C n×n , we can explicitly derive the following input to state transfer functions G n :

for n ≥ 1 and p + q = n. By multiplying with the output matrix C, we can further derive the input-output generalized frequency responses GFRFs (generalized transfer functions) as:

At this point, we are capable to write the Volterra series by using the above specific structure of the GFRFs that were derived with the growing exponential approach for the bilinear case. An important property to notice is that the n th kernel is a multivariate function of order n. It is obvious that the identification of the n th order involves an n-dimension frequency space. For that reason, next, we derive the general second symmetric kernel for the bilinear case with a 2-tone input. A double oscillatory input as:

where

In that case, with the growing exponential approach, we can write the input to state solution in steady state as:

6G p,q n = G( jω, ..., jω p-times ;jω, ...,jω q-times )

We are looking for the input to state frequency response G( jω 1 , jω 2 ). By substituting to the bilinear model in Eq. 19 and collecting the appropriate terms while at the same time using the symmetry G( jω 1 , jω 2 ) = G( jω 2 , jω 1 ), we conclude to:

where by using the resolvent notation and multiplying with C, we derive the 2-order generalized frequency response (symmetric kernel)

One way to measure Volterra kernels is by means of interpolation. This problem is equivalent to that of estimating the coefficients of a noisy polynomial. This interpolation scheme builds a linear system with a Vandermonde matrix. This matrix is invertible since the amplitudes are distinct and non-zero.

The inverse of a Vandermonde matrix can be explicitly computed and there are stable ways to solve these equations. The nth harmonic in the frequency domain is derived by applying a (single-sided) Fourier transform. More precisely, the explicit formulation is as follows:

We simplify the notation in order to reveal the adaptive method that will help us to estimate the GFRFs up a specific order. We can write the system which connects the harmonic information with the higher Volterra kernels as follows:

. . .

. . .

 By introducing the Hadamard product7 notation and simplify the δ's contribution to ones, we can compactly rewrite the above system in the following form:

The above system offers the level of approximation we want to achieve. Note that the frequency response Y depends on both the amplitude and the frequency, while the right hand side of Eq. 29 reveals the separation of the aforementioned quantities. As we neglect higher order Volterra kernels, the measurement set tends to be more and more corrupted by noise. Remark 4.1. Kernel separation and stage -approximation For a given system, the procedure consists in exciting it with an 1-tone input. By varying the driving frequency as well as, the amplitude, we can approximate GFRFs by minimizing the (2-norm) of the remaining over-determined systems.

7Hadamard product denotes as "•" and means that the multiplication is element wise.

where the n-"direction" gives us the threshold up to the specific harmonic that we would like to measure. The -"direction" give us the level of the kernel separation we want to achieve. For instance, for the 2nd stage approximation = 2 with α m ≈ 0, ∀m > = 2, and that results to n = 2 (up to 2nd harmonic).

The only difference between linear and bilinear models is the presence of the product between the input and the state that is scalled by the matrix N. As the (LF) is able to identify the linear part ( A,E,B,C) of the bilinear model the only thing that remains is the identification of the matrix N. The matrix N enters linearly in the following kernels:

• With a 1-tone input the H 1,1 2 can be written:

and the H 2,0 2 as:

• While with a 2-tone input the general H 2 can be written:

We introduce the following notations:

Note that, Eq. 33 be compactly rewritten as:

Assume that k measurements of the function H 2 are available (measured) for k different pairs (ω 1 , ω 2 ). By vectorizing in respect to the measurement set, we have:

Note that Eqs. (31,32,33) can be equivalently rewritten as one linear matrix equation, i.e., the one given in Eq. 36. By filling out the matrix [O ⊗ R] with the information from H 2 ( jω 1 ,jω 1 ) and from H 2 ( jω 1 , jω 1 ), has full rank. Hence, we are able to solve equation (36) and then identify the matrix N. All the symmetry properties of the kernels are appropriately used, e.g., conjugate-real symmetry. For n denoting the dimension of the bilinear model and k the number of measurements, we have the following two cases:

As the above lemma indicates, we need at least n 2 measurements to identify the matrix N from a bilinear system of dimension n.

To identify the n Volterra kernel, we need an n-tone input signal. As we want to identify the 2nd kernel, the input signal needs to be chosen as a double sinusoidal Eq. 24. The propagating harmonics are: e (j(m 1 -m 2 )ω 1 +j(m 3 -m 4 )ω 2 )t or more compactly e (±k jω 1 ±l jω 2 )t , where k, l ∈ N. The aim is to differentiate the (ω 1 + ω 2 ) harmonic from the others harmonics. By succeeding that, there will be no overlapping (commensurate frequencies) and we know what exactly we measure in the spectrum. More precisely, we want the following result to hold:

Suppose ω 2 = φω 1 , φ ∈ R. What are the suitable φ's where Eq. 38 holds?

By choosing φ that the equality in Eq. 39 doesn't hold, with harmonic mixing index as m = k + l (i.e., for 2nd stage approximation m = 2), it makes the harmonic (ω 1 + ω 2 ) uniquely defined in the frequency spectrum up to the m th harmonic.

To visualize this feature, we choose ω 1 = 1, and ω 2 = ω 1 φ = φ, for harmonic mixing index m = 2. Then, the constrains of φ in Fig. 4 are depicted with blue dots. Next, in Fig. 5 and on the left pane, one φ constraint that occurs commensurate harmonics is depicted, and on the right pane, one unique harmonic construction at (ω 1 + ω 2 ) is presented.

The next theoretical result allows us to construct sweeping frequency schemes to get enough measurements for the H 2 ( jω 1 , jω 2 ). So, for every ω 1 > 0 the following should hold :

where φ i are the constrains (see Fig. 4 blue dots).

With the framework that we have developed, we force the separation of ω 1 + ω 2 harmonic only under a specific mixing order m. We do not offer any general solution of the harmonic separation problem for multi-tone input, although techniques have been introduced such as in [13].

There, it was also stated that the solution of the full separation of harmonics is in general, not possible.

We start with a set of 1-tone inputs as u(t) = α 1 cos(ω i t), i = 1, ..., n, with α 1 < 1. For those n measurements we can estimate the linear kernel H 1 ( jω i ), the H 2 ( jω i , jω i ) and the H 2 ( jω i ,jω i ) by simply measuring the first harmonic as Y 1 , the second harmonic as Y 2 , and the DC term as Y 0 , from the frequency spectrum in Fig. 3. We can increase the accuracy of the method to -stage by varying the amplitude α as explained in section 4.2.

Since the (LF) reveals the underlying order of the system denoted with r, the value of n should be at least equal to 2r. Then, we can take the by analysing the singular value decay. Up to the previous step, we have identified the linear part with the (LF), and, we have filled the LS problem Eq. 36 with measurements from the diagonal of the 2nd kernel and from the the perpendicular to the diagonal axis (ω 1 , -ω 1 ). Those measurements contribute to the LS problem, but with an underdetermined (rank deficient) LS problem.

We need more measurements of H 2 to reach the full rank (r 2 ) solution that will lead to the identification of N. So, we proceed by measuring the H 2 off the diagonal with a 2-tone input as u(t) = α cos(ω 1,k t) + β cos(ω 2,k t), for a set of frequency pairs (ω 1 , ω 2 ) up to r 2 . The harmonic separation problem for the frequency (ω 1 , ω 2 ) appears now. To deal with this problem, we follow the solution proposed in section 4.4 (up to a mixing degree). Last, we solve the LS problem described in Eq. 36 by using all the symmetric properties of these kernels (i.e., real symmetry, conjugate symmetry and the fact that H 2 ( jω 1 , jω 2 ) = H 2 ( jω 2 , jω 1 )). An algorithm that summarizes the above procedure is presented below. Input: Use as control input the signals: u(t) = α cos(ω 1,k t)+ β cos(ω 2,k t), t ≥ 0, by sweeping the small amplitudes (< 1) and a particular range of frequencies.

Output: A reduced bilinear system of dimension-r: Σ r : (E r , A r , B r , C r , N r ) 1. Apply one-tone input u(t) with β = 0, ω 1,k for k = 1, . . . , n and collect the snapshots y(t) in steady state. Improve the measurements by solving system Eq. 29 by sweeping the amplitude for each frequency.

2. Apply Fourier transform for each frequency and collect the following k measurements:

• 2nd harmonic:

3. Apply the linear the (LF), see Algorithm 1 in [24] by using the measurements (e.g., H 1 ( jω 1,k ) = 2Y I ( jω 1,k )/α 1, for the 2nd stage approximation) and get the the order r reduced linear model. If the second harmonic and higher harmonics are equal with zero, the underlying system is linear so N is the zero matrix 0 r×r as well. Otherwise the system is nonlinear and computing the bilinear matrix N will improve the accuracy.

4. Apply the 2-tone input u(t) = α cos(ω 1,k t) + β cos(ω 2,k t) to get enough measurements (≤ r 2 ) to produce a full rank LS problem. Measure the (ω 1 + ω 2 ) harmonic as explained in section 4.4 and get the estimations for the 2nd kernel as:

5. Solve the full rank least square problem described in Eq. 36 and compute the real-valued bilinear matrix N.

Example 5.1 (A small bilinear identification). The aim in this example is to identify the bilinear model from time domain measurements. Consider the following controllable/observable bilinear model Eq. 19 of dimension-2 with a non-symmetric matrix N, zero initial condition and matrices as:

We simulate the system in time domain with a single sinusoidal input u(t) = 0.01 cos(2πωt), by varying the frequency ω ∈ [0.5 1 1.

With the estimations of the linear transfer function and by using the (LF) as the data-driven identification and reduction tool for linear systems, we identify the linear system ( Ẽ, Ã, B, C). We stopped at the 4th measurement due to the fact that the underlying system is of second order (McMillan degree 2). Otherwise, more measurements will be needed to have a reviling decay of the singular values as in the next Fig. 6. The singular values decay offers a choice of reduction. As long as the simulation of the system is done, with time step dt = 1e -4, the singular values with magnitude below that threshold are neglected. Construction of the linear system with order r = 2 by using the theoretical noise-free measurements (subscript "t") appears next: This example illustrating the bilinear modeling and reduction concepts proposed in [5] for the viscous Burger's equation from time domain simulations. We kept the same set up, as in the introductory example and we present the corresponding results with initial system dimension n = 110 reduced by the proposed method to order r = 2 with the first neglected singular value to be σ 3 /σ 1 = 2.576 • 10 -5 . In Fig. 9, evaluation results are presented. Last, in Fig. 10,

What makes this algorithm feasible is the combination of the data-driven Loewner framework with the non-linear Volterra series framework. The proposed method offers approximate bilinear system identification from time domain measurements, since it is not possible to measure the corresponding kernels exactly. Our proposed method uses only input-output measurements without requiring statespace access. The evaluations of the kernels and the outputs (y l : linear, ỹb reduced bilinear (r=2)) took place over the domains that depicted in Figs. (8,7,9,10).