Mathematical models are commonly used to simulate, optimize, and control the behavior of real dynamical processes. A common way to derive those models is to use the first principles, generally leading to a set of ordinary or partial differential equations. For high complex dynamics, fine discretization leads to high fidelity models, which require numerous equations and variables. In some situations, the high model is given as a black box setup, i.e., by solvers that allow the computation of the full model states for a given set of initial conditions and inputs, but does not provide the dynamical system realization. In order to better understand such dynamics processes, it is often beneficial to construct surrogate models using simulated data. This justifies the development of identification or data-driven model reduction methods. Indeed, with the ever-increasing availability of measured/simulated data in different scientific disciplines, the need for incorporating this information in the identification and reduction process has steadily grown. The data-driven model reduction problem consists of determining low-order models from the provided data obtained either by experimentation or numerical simulations. Methods such as Dynamic Mode Decomposition (DMD) have drawn considerable research endeavors.
DMD is a data-driven method for analyzing complex systems. The purpose is to learn/ extract the important dynamic characteristics (such as unstable growth modes, resonance, and spectral properties) of the underlying dynamical system by means of measured time-domain data. These can be acquired through experiments in a practical setup or artificially, through numerical simulations (by exciting the system). It was initially proposed in  in the context of analyzing numerical and experimental fluid mechanics problems. Additionally, it is intrinsically related to the Koopman operator analysis, see [26, 21]. Since its introduction, several extensions have been proposed in the literature, e.g., the exact DMD , the extended DMD  and the higher-order DMD . Also, in order to address control problems, DMD with control inputs was proposed in , and then extended to the case where outputs are also considered in [1, 9]. The reader is refereed to  for a comprehensive monograph on the topic.
Often, nonlinear systems modeling physical phenomena have a particular known structure, such as bilinear and quadratic terms. In the present work, our primary goal is to embed nonlinear structures in the DMD framework. To this aim, we propose an identification and data-driven reduction method based on the classical DMD approach allowing to fit a bilinear and quadratic-bilinear structures to the measured data. The choice to fit such terms is due to the fact most systems with analytical nonlinearities (e.g., rational, trigonometrical, polynomial) can be exactly reformulated as quadratic-bilinear systems . Our work is rooted in the two variants, DMD with control and input-output DMD, and can be considered as an extension of those methodologies.
There exist vast literature on learning nonlinear dynamics from data, and we review the most relevant literature for our work. One approach is the so-called Loewner framework, which enables to construct low-order models from frequency domain data. It was initially proposed in, and later extended to bilinear  and quadratic-bilinear case . Another approach is the operator inference, proposed in . This approach infers polynomial low-order models as a solution of a least-squares problem based on the initial conditions, inputs, trajectories of the states. This approach was recently extended to systems with non-polynomials . Also, the authors in  show how the use of lifting transformations can be beneficial to identify the system. Finally, the approach proposed in , introduces a method based on operator inference enabling to learn exactly the reduced models that are traditionally constructed with model reduction. However, to the best of our knowledge, there exists no available methodology incorporating bilinear and quadratic-bilinear structures to DMD.
The rest of the paper is organized as follows. Section 2 recalls some results on the classical DMD, DMD with control and the input-output DMD. In Section 3 we present the main contribution of the paper, which is the incorporation of bilinear and quadratic-bilinear terms in the DMD setup. Finally, in Section 4, we demonstrate the proposed methodology for different examples, such as the Burgers’ equation and the coupled van der Pol oscillators.
2 Dynamic mode decomposition
In this section, we briefly recall the classical DMD framework 
. To this aim, we analyze time-invariant systems of ordinary differential equations (ODEs) written compactly in a continuous-time setting as follows
is the state vector andis the system nonlinearity.
By means of sampling the variable in (1) at uniform intervals of time, we collect a series of vectors for sampling times . For simplicity, denote .
DMD aims at analyzing the relationship between pairs of measurements from a dynamical system. The measurements and , as previously introduced, are assumed to be approximately related by a linear operator for all .
where . This approximation is assumed to hold for all pairs of measurements. Next, group together the sequence of collected snapshots of the discretized state and use the following notations:
The DMD method is based on finding a best-fit solution of an operator A so that the following relation is (approximately) satisfied
which represents the block-version of equation (2). Moreover, the above relation does not need to hold exactly. Previous work has theoretically justified using this approximating operator on data generated by nonlinear dynamical systems. For more details, see . A best-fit solution is explicitly given as follows:
where is the Moore-Penrose pseudo-inverse of matrix . By best fit we mean the solution that minimizes the least-squares error in the Frobenius norm (see ).
The so-called DMD modes are given by the eigenvectors of matrixin (5), collected in matrix with . These spatial modes of system (1) are computed at a single frequency and are connected to the Koopman operator, see .
In this work, we will mainly focus on the construction of the reduced-order model rather than the evaluation of the DMD modes.
2.1 Dynamic mode decomposition with control (DMDc)
Dynamic mode decomposition with control (DMDc) was introduced in  and it modifies the basic framework characterizing DMD. The novelty is given by including measurements of a control input . It is hence assumed that the dynamics of the original system of ODEs includes an input dependence, i.e.
In this setup, a trio of measurements are now assumed to be connected. The goal of DMDc is to analyze the relationship between a future state measurement with the current measurement and the current control .
The motivation for this method is that, understanding the dynamic characteristics of systems that have both internal dynamics and applied external control is of great use for many applications, such as for controller design and sensor placement.
The DMDc method is used to discover the underlying dynamics subject to a driving control input by quantifying its effect to the time-domain measurements corresponding to the underlying dynamical system.
A pair of linear operators represented by matrices and provides the following dependence for each trio of measurement data snapshots
Next, denote the sequence of control input snapshots with
The first step is to augment the matrix with the row vector and similarly group together the and matrices by using the notations:
The matrix introduced above will be referred to as the system matrix since it incorporates the matrices corresponding to the system to be fitted.
By letting the index k vary in the range , one can compactly rewrite the equations in the following matrix format:
Thus, similar to standard DMD, compute a pseudo-inverse and solve for as
To solve (11
), we compute a singular value decomposition (SVD) of the augmented data matrixas follows
where the full scale and reduced-order matrices have the following dimensions
The truncation index is denoted with , where . The pseudo-inverse is computed using the matrices from the SVD in (12), i. e., as .
By splitting up the matrix as , recover the system matrices as
As mentioned in, there is one additional step. By performing another (short) SVD of the matrix , write
where . Note that the two SVDs will likely have different truncation values. The following reduced-order approximations of and are hence computed as:
2.2 Input-output dynamic mode decomposition
In this section we discuss the technique proposed in  known as input-output dynamic mode decomposition (ioDMD). This method constructs an input-output reduced-order model and can be viewed an extension of DMDc for the case with observed outputs. As stated in the original work , this method represents a combination of POD and system identification techniques. The proposed method discussed here is similar in a sense to the algorithms for subspace state space system identification (N4SID) introduced in  and can be also applied to large-scale systems.
We consider as given a system of ODEs whose dynamics is described by the same equations as in (6). Additionally, assume that observations are collected in the variable , as function of the state variable and of the control , written as
As before, the next step is to collect snapshots of both variable and of the output sampled at some positive time instances . Again, for simplicity of the exposition, denote with .
We enforce the following dependence for each trio of measurement data snapshots given by
Afterwards, collect the output values in a row vector as follows
The ioDMD method aims at fitting the given set of snapshot measurements collected in matrices and vectors and to a linear discrete-time system characterized by the following equations
where as before and , and also . Note that the first equation in (19) exactly corresponds to the driving matrix equation of DMDc presented in (11). Moreover, write the system of equations in (19) compactly as
Next, we adapt the definition of the system matrix from (9) by incorporating an extra line as follows
while is as before. Introduce a new notation that will become useful also in the next sections. It represents an augmentation of the shifted state matrix with the output observations vector , i.e.,
Again, the solution of equation (20) will be computed as a best-fit type of approach. Hence, similarly to the DMDc case, recover the matrices , , and by computing the pseudo-inverse of matrix and writing
3 The proposed extensions
In this section, we present the main contribution of the paper. We propose extensions of the methods previously introduced in Sections 2.1 and 2.2, e.g. DMDc and, respectively ioDMD to fit nonlinear structured systems. More specifically, the discrete-time models that are fitted using these procedures will no longer be linear as in (19); the new models will contain nonlinear (bilinear or quadratic) terms.
3.1 Bilinear Systems
Bilinear systems are a class of mildly nonlinear systems for which the nonlinearity is given by the product between the state variable and the control input. More exactly, the characterizing system of ODEs is written as in (6) but for a specific choice of mapping , i.e. . Additionally, assume that the observed output depends linearly on the state . Hence, it what follows, we will make use of the following description of bilinear systems (with a single input and a single output)
where the matrix scales the product of the state variable with the control input .
In practice, bilinear control systems are used to approximate nonlinear systems with more general, analytic nonlinearities. This procedure is known as Carleman’s linearization; for more details see .
Bilinear systems are a class of nonlinear systems that received considerable attention in the last four or five decades. Contributions that range from realization theory in , classical system identification in  or to subspace identification in . In more recent years (last two decades), model order reduction of bilinear systems (in both continuous and discrete-time domain) was extensively studied with contributions covering balanced truncation in , Krylov subspace methods in 
, interpolation-basedmethod in [4, 5], or data-driven Loewner approach in .
3.1.1 The general procedure
We start by collecting snapshots of the state for multiple time instances . We enforce that the snapshot at time depends on the snapshot in the following way
where . One can hence write , with and then introduce the matrix as
The next step is to augment the matrix with matrix and denote this new matrix with as an extension of the matrix previously introduced in (9), i.e.,
For the case in which we extend the DMDc method in Section 2.1 to fitting bilinear dynamics (no output observations), we propose a slightly different definition for the matrix . We hence append the matrix to the originally introduced system matrix in (9). Then, equation (26) can be written in a factorized way as , where the matrices for this particular setup are as follows
Alternatively, for the case where output observations are also available, we enforce a special bilinear dependence for each trio of measurement data snapshots as
Afterwards, we collect the equations in (30) for each index and hence write
This system of equations can be then written in a factorized way as before, i.e., , where the matrices for this particular setup are given below
Finally, the last step is to recover the matrix and split it block-wise in order to put together a system realization. Consequently, this all boils down to solving the equation (in either of the two cases, with or without output observations included).
As shown in the previous sections, e.g., in (11), solving for is based on computing the pseudo-inverse of matrix . More precisely, we write the solution as
3.1.2 Computation of the reduced-order matrices
In this section we present specific/practical details for retrieving the system matrices in the case of the proposed procedure in Section 3.1.1. We solve the equation for which the matrices are given as in (33), i.e. the case containing output observations. We compute an SVD of the augmented data matrix giving
where the full-scale and reduced-scale matrices derived from SVD are as follows
The truncation index is denoted with , where .The computation of the pseudo-inverse is done by the SVD approach, i.e., By splitting up the matrix as , one can recover the system matrices as
By performing another (short) SVD for the matrix , we can write
where . Note that the two SVDs could have different truncation values denoted with p and r. Using the transformation , the following reduced-order matrices can be computed:
3.1.3 Conversions between discrete-time and continuous-time representations
The DMD-type approaches available in the literature identify continuous-time systems by means of linear discrete-time models. In this contribution, we make use of the same philosophy, in the sense that the models fitted are discrete-time. We extend the DMDc and ioDMD approaches by allowing bilinear or quadratic terms to appear in the these models as well.
As also mentioned in , one can compute a continuous-time model that represents a first-order approximation of the discrete time model obtained by DMD-type approaches.
Assume that we are in the bilinear setting presented in Section 3.1 and that we already have computed a reduced-order discrete-time model given by matrices , i.e., following the explicit derivations in (38). Then, a continuous-time model can also be derived. By assuming that the standard first-order Euler method was used for simulating the original system (with a small enough time step-size ), we can write that
Observe that for the ioDMD-type of approaches, the feed-through terms that appear in the output-state equation are the same in both discrete and continuous representations.
3.2 Quadratic-Bilinear Systems
Next, we extend the method in Section 2.1 for fitting another class of nonlinear systems, i.e. quadratic-bilinear (QB) systems. Additional to the bilinear terms that enter the differential equations, we assume that quadratic terms are also present. More precisely, the system of the ODEs is written as in (6) but for a specific choice of nonlinear mapping , i.e.,
where ’’ denotes the Kronecker product, the matrix scales the product of the state with itself, and is as in Section 3.1.
Quadratic-bilinear systems appear in many applications for which the original system of ODEs inherently has the required quadratic structure. For example, after semi-discretizing the Burgers’ or Navier-Stokes equations in the spatial domain, one obtains a system of differential equations with quadratic nonlinearities (and also with bilinear terms). Moreover, many smooth analytic nonlinear systems, that contain combinations of nonlinearities such as, e.g. exponential, trigonometric, polynomial functions, etc. can be equivalently rewritten as QB systems. This is performed by employing so-called lifting techniques. More exactly, one needs to introduce new state variables in order to simplify the nonlinearities and hence derive new differential equations corresponding to these variables. Model order reduction of QB systems was a topic of interest in the last years with contributions ranging from projection-based approaches in [15, 6] to optimal -based approximation in , or data-driven approaches in the Loewner framework in [14, 2].
Similarly to the procedure described in Section 3.1, we enforce that the snapshot at time depends on the snapshot in the following way
Next, by varying the in the range , compactly rewrite the equations in (40) in the following matrix format:
with and Here, is the unit vector of length that contains the 1 on position . Additionally, we introduce the matrix that depends on the state matrix as
Note that the equality holds as follows . Next, we augment the matrix with both matrices and and group together the matrices and by using the notations:
Hence, by using these notations, it follows that we can rewrite the equation (41) as follows
Thus, similarly to the previous methods covered so far, we can recover the matrix , and hence system matrices and ), by means of a least-squares approach that involves computing pseudo-inverse matrix, i.e., .
As previously shown in Section 3.1, we can again adapt this proposed procedure for fitting QB systems in the ioDMD format by involving output observations measurements . The procedure for quadratic-bilinear systems is similar to that for bilinear systems and we prefer to skip the exact description to avoid duplication. For more details, see the derivation in Section 6.1.
4 Numerical Experiments
4.1 The viscous Burgers’ equation
Consider the partial differential viscous Burgers equation:
with i.c. . The viscosity parameter is denoted with .
The Burgers’ equation has a convective term, an unsteady term and a viscous term; it can be viewed as a simplification of the Navier-Stokes equations.
By means of semi-discretization in the space domain, one can obtain the following nonlinear (quadratic) model (see ) described by the following system of ODEs
Next, by means of the Carleman linearization procedure in , one can approximate the above nonlinear system of order with a bilinear system of order . The procedure is as follows: let be the original state variable in (43). Then, introduce the augmented state variable corresponding to the system described by the following equations
The continuous-time bilinear model in (44) is going to be used in following the numerical experiments.
Start by choosing the viscosity parameter to be . Then, choose as the dimension of the original discretization and hence, the bilinear system in (44) is of order . Perform a time-domain simulation of this system by approximating the derivative as follows . We use as time step and the time horizon to be . The control input is chosen to be .
Hence, collect snapshots of the trio that are arranged in the required matrix format as presented in the previous sections. The first step is to perform an SVD for the matrix . The first 200 normalized singular values are presented in Figure 1. Choose the tolerance value which corresponds to a truncation order of (for computing the pseudo-inverse of matrix ). On the same plot in Figure 1 we also display the normalized singular values of matrix . Note that machine precision is reached at the singular value. We select three tolerance values for truncating matrices obtained from the SVD of .
Start by choosing the first tolerance value, e.g. . This corresponds to a truncation value of . We compute term and a also bilinear feed-through term with . We simulate both the original large-scale bilinear system and the reduced-order system. The results are presented in Figure 2. Note that the observed output curves deviate substantially. We hence need to decrease the tolerance value.
For the next experiment, choose the tolerance value to be . This corresponds to a truncation value of . After computing the required matrices, notice that the D term is again numerically 0, while the norm of the matrix slightly decreases to the value . Perform numerical simulations and depict the two outputs and the approximation error in Figure 3. Observe that the approximation quality significantly improved, but there is still room for improvement.
Finally, the tolerance value is chosen as . For this particular choice, it follows that the truncation value is . In this case, the output of the reduced-order model faithfully reproduces the original output, as it can be observed in Figure 4. Note also that the approximation error stabilizes within the range .
4.2 Coupled van der Pol oscillators
Consider the coupled van der Pol oscillators along a limit cycle example given in . The dynamics are characterized by the following six differential equations with linear and nonlinear (cubic) terms:
Choose the output to be . Hence the state-output equation is written as with and . Choose the parameters in (45) as follows: and .
Note that by introducing three additional surrogate states, e.g. and , one can rewrite the cubic nonlinear system in (45) of order as an order quadratic-bilinear system.
Perform time-domain simulations of the cubic system of order and collect data from 500 snapshots using the explicit Euler method with step size . The chosen time horizon is hence [0,5]s. The control input is a square-wave with amplitude 30, i.e., .
Compute the pseudo-inverse of matrix and select as truncation value (the 20th normalized singular value drops below machine precision).
We compute a reduced-order quadratic-bilinear model of order . We made this choice since the fifth normalized singular value of matrix is 5.8651e-04 while the sixth is numerically 0, i.e., 3.5574e-16. We hence fit an order quadratic-bilinear system that approximates the original order cubic polynomial system. Note that the only non-zero feed-through quantity in the recovered state-output equation si given by .
The time-domain simulations are depicted in Fig. 5. One can observe that the two outputs match quite well. In this particular setup, it follows that the response of the 6th order cubic system (that can be equivalently written as a 9th order QB system) can be accurately approximated with the response of a 5th order QB system. The approximation error is also presented in Fig. 5.
In this paper, we have proposed extensions of the DMDc and ioDMD recently proposed methods. The philosophy is similar to that of the original methods, but instead of fitting discrete-time linear systems, we impose a more complex structure to the fitted models. More precisely, we fit bilinear or quadratic terms to augment the existing linear quantities (both in the differential and in the output equations). The numerical results presented were promising, and they have shown the strength of the method. Indeed, there is a clear trade-off to be made between approximation quality and the dimension of the fitted model.
Nevertheless, this represents a first step towards extending DMD-type methods, and a more involved analysis of the method’s advantages and disadvantages could represent an appealing future endeavor. Moreover, another contribution could be made by comparing the proposed methods in this work with the recently-introduced operator inference-type methods. For the quadratic-bilinear case, additional challenges arise when storing the large-scale matrices involved and also when computing the classical SVD for such big non-sparse matrices.
6.1 Computation of the reduced-order matrices for the quadratic-bilinear case
In this section we present practical details for retrieving the system matrices in the case of the proposed procedure in Section 3.2. We solve the equation for which the matrices are given as in (42), i.e. the case without output observations. We again utilize an SVD, now performed on the matrix , i.e.,
where the full scale and reduced SVD matrices have the following dimensions
The truncation index is denoted with , and write as before
By splitting up the matrix as , with
recover the matrices
Again, perform an additional SVD, e.g. , where . Using the transformation , the following reduced-order approximations are computed:
-  (2016) Wind farm flow modeling using an input-output reduced-order model. In 2016 American Control Conference (ACC), pp. 506–512. Cited by: §1, §2.2.
-  (2018) On the Loewner framework for model reduction of Burgers’ equation. In Active Flow and Combustion Control, R. King (Ed.), pp. 255–270. Cited by: §3.2.
-  (2016) Model reduction of bilinear systems in the Loewner framework. SIAM Jour. Sci. Comp. 38 (5), pp. B889–B916. Cited by: §1, §3.1.
-  (2011) Generalised tangential interpolation for model reduction of discrete-time mimo bilinear systems. International Journal of Control 84 (8), pp. 1398–1407. Cited by: §3.1.
-  (2012) Interpolation-based -model reduction of bilinear control systems. SIAM Journal on Matrix Analysis and Applications 33 (3), pp. 859–885. Cited by: §3.1.
-  (2015) Two-sided projection methods for nonlinear model order reduction. SIAM Journal on Scientific Computing 37 (2), pp. B239–B260. Cited by: §3.2, §4.1.
-  (2018) -Quasi-optimal model order reduction for quadratic-bilinear control systems. SIAM Journal on Matrix Analysis and Applications 39 (2), pp. 983–1032. Cited by: §3.2.
-  (2020) Operator inference for non-intrusive model reduction of systems with non-polynomial nonlinear terms. arXiv preprint arXiv:2002.09726. Cited by: §1.
-  (2018) On reduced input-output dynamic mode decomposition. Adv. Comput. Math 44, pp. 1751–1768. Cited by: §1, §2, §3.1.3.
-  (2010) Krylov subspace methods for model order reduction of bilinear control systems. Systems Control Letters 59 (8), pp. 443–450. Cited by: §3.1.
-  (2012) Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses. Journal of nonlinear science 22 (6), pp. 887–915. Cited by: §1.
-  (1990) A method for bilinear system identification. Vol. IFAC Proceedings23, pp. 143–148. Cited by: §3.1.
-  (1999) Subspace identification of bilinear systems subject to white inputs. IEE Trans. Autom. Ctrl. 44 (6), pp. 1157–1165. Cited by: §3.1.
-  (2018) Data-driven model order reduction of quadratic-bilinear systems. Numerical Linear Algebra with Applications 25 (6), pp. e2200. Cited by: §1, §3.2.
-  (2011) A projection-based nonlinear model order reduction approach using quadratic-linear representation of nonlinear systems. IEEE Trans. Comput.-Aided Design Integr. Circuits Syst 30, pp. 1307–1320. Cited by: §1, §3.2.
-  (1973) Realization theory of bilinear systems.. In Geometric Methods in System Theory, D. Q. Mayne and R. W. Brockett (Eds.), Cited by: §3.1.
-  (2017) Empirical differential balancing for nonlinear systems. Vol. 20th IFAC World Congress, Toulouse, France, July 9-1450, pp. 6326–6331. Cited by: §4.2.
-  (2016) Dynamic mode decomposition: data-driven modeling of complex systems. SIAM. Cited by: §1.
-  (2017) Higher order dynamic mode decomposition. SIAM Journal on Applied Dynamical Systems 16 (2), pp. 882–925. Cited by: §1.
-  (2007) A framework for the solution of the generalized realization problem. Linear Algebra and its Applications 425 (2-3), pp. 634–662. Cited by: §1.
-  (2013) Analysis of fluid flows via spectral properties of the koopman operator. Annual Review of Fluid Mechanics 45, pp. 357–378. Cited by: §1, §2.
-  (2016) Data-driven operator inference for nonintrusive projection-based model reduction. Computer Methods in Applied Mechanics and Engineering 306, pp. 196–215. Cited by: §1.
-  (2019) Sampling low-dimensional markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference. arXiv preprint arXiv:1908.11233. Cited by: §1.
-  (2016) Dynamic mode decomposition with control. SIAM Journal Applied Dynamical Systems 15 (1), pp. 142–161. Cited by: §1, §2.1.
-  (2019) Transform & learn: a data-driven approach to nonlinear model reduction. In AIAA Aviation 2019 Forum, pp. 3707. Cited by: §1.
-  (2009) Spectral analysis of nonlinear flows. Journal of Fluid Mechanics 641, pp. 115–127. Cited by: §1.
-  (1981) Nonlinear system theory - the volterra/wiener approach. University Press. Cited by: §3.1, §4.1.
-  (2010) Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, pp. 5–28. Cited by: §1, §2.
-  (2014) On dynamic mode decomposition: Theory and applications. J. Comput. Dynam. 1, pp. 391–421. Cited by: §2.
-  (2014) On dynamic mode decomposition: theory and applications. Journal of Computational Dynamics 1, pp. 391–421. Cited by: §1.
-  (1994) N4SID: subspace algorithms for the identification of combined deterministic-stochastic systems. Automatica 30 (1), pp. 75–93. Cited by: §2.2.
-  (1999) On gramians and balanced truncation of discrete-time bilinear systems. International Journal of Control 76 (2003), pp. 414–427. Cited by: §3.1.