Krylov Subspace Methods for Model Reduction of Quadratic-Bilinear Systems

The authors propose a two sided moment matching method for model reduction of quadratic-bilinear descriptor systems. The goal is to approximate some of the generalised transfer functions that appear in the input–output representation of the non-linear system. Existing techniques achieve this by utilising moment matching for the first two generalised transfer functions. In this study, they derive an equivalent representation that simplifies the structure of the generalised transfer functions. This allows them to extend the idea of two sided moment matching to higher subsystems which was difficult in the previous approaches. Numerical results are given for some benchmark examples of quadratic-bilinear systems.


Introduction
Consider a multi-input multi-output (MIMO) quadratic-bilinear descriptor system of the form Here , , ∈ ℝ × , ∈ ℝ × 2 , ∈ ℝ × and ∈ ℝ × are the state-space matrices and () ∈ ℝ , () ∈ ℝ and () ∈ ℝ are the state, input and output vectors, respectively.This class of nonlinear systems has applications in many areas including the simulation of fluid flow, electrical circuits and some biological systems, cf.[1].Also, a large class of non-linear systems can be written as quadratic-bilinear differential algebraic equations (QBDAEs) by utilising exact transformations [2], extending its application areas even further.
In this paper, we discuss the approximation of QBDAEs by constructing reduced order models.That is, we construct with ^() ∈ ℝ and ≪ such that ^() is close to () in an appropriate norm.Recently, projection based (generalised-) moment matching techniques [2][3][4] have been used in the literature to construct such reduced order systems, in contrast to trajectorybased methods such as proper orthogonal decomposition (POD) [5], the reduced basis method, POD with discrete empirical interpolation method [6], and, the trajectory piecewise linear method [7].Since all these trajectory based methods share the disadvantage of input dependency, generalised moment matching techniques are particularly useful for systems with varied input function, which is common in control and optimisation problems.
For details on non-linear model reduction techniques, we refer a recent survey paper [8].
Projection involves identifying suitable basis matrices and (the columns of each matrix span a particular subspace), approximating () by ^() and ensuring Petrov-Galerkin conditions.This leads to the following reduced state-space matrices: ^= T , ^= T , ^= T ( ⊗ ), ^ = T , = 1, …, , ^= T , ^= . ( Clearly the reduced system matrices depend on the choice of the basis matrices and for a given Σ.In moment matching techniques, these projection matrices are constructed such that the reduced system matches some of the generalised moments associated with the underlying generalised transfer functions of the QBDAE system at fixed interpolation points.In [2], the procedure of moment matching is restricted to orthogonal projection (i.e.= ) with the moments matched for each of the first two generalised transfer functions.This idea has been extended recently to oblique projection [4] which improves the approximation quality of the reduced system, but because of the complex structure of the moments corresponding to third and higher dimensional generalised transfer functions, the higher dimensional moments are again ignored.
In this paper, we utilise the known connection between different forms of the generalised transfer functions in order to identify an equivalent representation of the generalised transfer functions which simplifies the structure of the moments.This allows us to use two sided projection based moment matching techniques also for higher dimensional transfer functions.
In Section 2, we briefly overview some background theory and existing moment matching techniques for model reduction of QBDAEs.To extend the idea of moment matching to higher subsystems, we derive a simplified form of multivariate transfer functions in Section 3. Based on the simplified form, we propose our moment matching framework using orthogonal as well as oblique projections.These results are given in Section 4. Numerical results are presented in Section 5 and, finally, the conclusions and future work are given in Section 6.

Background
The input-output representation for single input quadratic-bilinear systems can be expressed by the Volterra series expansion of the output () with quantities analogous to the standard convolution operator.That is (see ( 4)) where it is assumed that the input signal is one-sided, () = 0 for < 0. In addition, each of the generalised impulse responses, ℎ ( 1 , …, ), also called thedimensional kernel of the subsystem, is assumed to be one-sided.In terms of the multivariable Laplace transform, the -dimensional subsystem can be represented as where ( 1 , …, ) is the multivariable transfer function of thedimensional subsystem.The above equation follows by applying the convolution property of the multivariable Laplace transform to (4), see [9] for details.If the multivariable Laplace transforms of all subsystems, that are ( 1 , …, )'s, and the input, (), are known, then the inverse Laplace transforms can be computed to identify ( 1 , …, ).The output () becomes The generalised transfer functions in the output expression (5) are in what is called the triangular form [9] of the functions.We denote the -dimensional triangular form by tri [] ( 1 , …, ).There are some other useful forms such as the symmetric and the regular forms of the multivariable transfer functions as discussed in [9].The triangular form is related to the symmetric form by the following expression: where the summation includes all !permutations of 1 , …, .Also, the triangular form can be connected to the regular form of the transfer function by using According to [9], the structure of the generalised symmetric transfer functions can be identified by the growing exponential approach.The structure of these symmetric transfer functions for the first three subsystems of the quadratic-bilinear system (1) can be written as (see (9)) .Clearly, if a model reduction approach can ensure that with ^sym [] ( 1 , …, ) being the multivariate transfer functions of the reduced system Σ ^, we can expect that the output () is well approximated by ^().This idea was initially utilised in [2] to construct a reduced quadratic-bilinear system with orthogonal projection of the first two subsystems using the symmetric transfer functions.Recently, this approach was extended in [4] to the oblique projection framework in order to improve the quality of the reduced model.The complex structure of the third and higher symmetric transfer functions has again restricted these projection techniques to the moment matching of the first two subsystems only.In the following we briefly review the oblique projection framework.Before proceeding further, we discuss some properties and notations: where , , ∈ ℝ are arbitrary and it is assumed that ( ⊗ ) = ( ⊗ ) holds.The matrix (2) is the 2- matricisation of the three-dimensional tensor ℋ ∈ ℝ × × having as its 1-matricisation, see [10].To recycle vectors for approximation subspaces, it is assumed in [4] With these settings, the second symmetric transfer function becomes sym [2] The following summarises the result introduced in [4].Lemma 1: Let ∈ ℂ be the interpolation points and ∉ {Λ(, ), Λ( , )}, where Λ(, ) represents the generalised eigenvalues of the matrix pencil − .Assume that ^= T is non-singular and ^, ^, ^, ^, ^ are as in (3) with full rank matrices , ∈ ℝ × such that (see equation below).Then the reduced QBDAE satisfies the following (Hermite) interpolation conditions: (see equation below).Proof: See [4] for a proof.□ In the remaining part of this paper, our goal is to identify the regular form of the multivariate transfer functions that can hopefully simplify the moment matching concept and allows us to use this new framework for third and higher subsystems.

Regular form of multivariate transfer functions
In this section, we utilise the connections between different forms of the multivariate transfer functions discussed in the previous section to identify the regular form of the corresponding functions.
Using (7), we observe that the equivalent symmetric form is the same as in (9).Similarly for higher values of , one can show that the regular form in (13) holds.□ Remark 1: The regular and triangular forms include − 1 sums of Kronecker products which is much smaller as compared to the corresponding symmetric form.Also in the symmetric form, it is difficult, if not impossible, to represent a general th-dimensional multivariate transfer function.
Remark 2: The symmetric form is exactly equal to the triangular form if we assume that 1 = ⋯ = = and the two forms are equal to the regular form if the regular variables are 1 = , 2 = 2, …, = .

Multimoment-matching with the regular form
In this section, we propose orthogonal as well as oblique projection techniques for model reduction of quadratic-bilinear systems using the regular form.We begin with the case of multimoment-matching for the first two subsystems only.The general case, where the multimoments associated with third and higher subsystems are also matched, is discussed later.
The above transfer functions are similar to symmetric transfer functions, if in the symmetric case 1 = 2 = , which is assumed in Lemma 1.Thus, the framework for interpolation of the transfer functions in the regular form is similar to Lemma 1.However, the interpolation of the partial derivatives with respect to 1 or 2 varies.The following theorem shows the interpolation conditions in the regular case for the first two subsystems.Theorem 2: Let ∈ ℂ be the interpolation points and × ∉ {Λ(, ), Λ( , )}, = 1, …, .Assume that ^= T is non-singular and ^, ^, ^, ^, ^ are as defined in (3) with , ∈ ℝ × such that (see ( 15)) Then the reduced QBDAE ensures that the following holds: reg [1] ( ) = ^reg [1] ( ), reg [2] ( , 2 ) = ^reg [2] ( , 2 ) and, in addition, interpolates all combinations of the multimoments that can be written as T and T ( + ( ⊗ )), where ∈ span() and ∈ span().Proof: Equation (16) holds due to the structure of and its proof is as in Lemma 1.To prove the second part, let ^∈ ℝ × and ^∈ ℝ × be defined such that (see (17)) Also let ^∈ span( ^) and ^∈ span( ^).Assume that , ^ (and , ^ ) are the same linear combination of the columns of the matrices on the right-hand side of ( 15) and (17), respectively.Then analogous to the discussion of one-sided projection, it is easy to show that This means that Similarly, T ( + ( ⊗ )) = ^T( ^ ^+ ^( ^⊗ ^)).□ The complete approach of one-sided projection for interpolatory model reduction of quadratic-bilinear systems, using regular generalised transfer functions, is shown in algorithm 1 (Fig. 1).

General case
The transfer function of the th subsystem in the regular form with the assumption that 1 = , 2 = 2, …, = can be written as The transfer function shows that one can recycle vectors in the construction of the projection matrix.This is shown in the following lemma which extends the orthogonal projection technique [2] for mutimoment matching to third and higher subsystems.

Numerical results
We demonstrate our results with two benchmark examples.These include the non-linear RC circuit example [1] and the 1D Burgers' equation [11].In each case, we use the iterative rational Krylov algorithm [12] on the corresponding linear part of the example to identify Süditaliens interpolation points.

Non-linear RC circuit
We consider a non-linear RC circuit as shown in Fig. 3, which is a benchmark example for non-linear model reduction [2,4] The non-linearity in the system is due to the diode I-V characteristics, given by () = 40 − 1, where is the node voltage.The current is treated as the input and the voltage 1 () at node 1 as the output of the system.Using Kirchhoff's current law at each of the nodes and assuming a normalised capacitance, = 1, we have where (()) is the non-linear function and = T is the first column of the × identity matrix.This non-linear model can be transformed [2] to an equivalent quadratic-bilinear descriptor system with size = 2.We choose = 500, which results in a quadratic-bilinear system of order 1000.
We reduce the order of the QBDAE system to = 10 by using orthogonal as well as oblique projections such that interpolation for the first two subsystems with the regular form of the transfer functions is ensured.The results are compared with the existing interpolation results based on the symmetric form, for the exponential function − as system input.The output response and the relative error are shown in Fig. 4a and b, respectively.
As discussed in Remark 2, the symmetric form is equal to the regular form under some interpolation conditions.Since we are using these interpolation conditions in Section 4, the projection matrix will not change in the symmetric and the regular forms.Thus the one sided projection results are the same in both the symmetric and regular forms.The partial derivatives of the regular and symmetric forms with respect to 1 and 2 are, however, different and therefore the oblique projection matrix is different in the symmetric and regular forms.As shown in Figs.4a and b, the relative error associated with the proposed two-sided projection technique in the regular case is comparable to the existing twosided projection technique in the symmetric case.
Interpolation of third and higher subsystems requires some simplifications in order to use the two-sided symmetric projection approach.Since the primary source of two-sided symmetric projection approach [4] is restricted to the first two subsystems, for interpolation of higher subsystems, we are not comparing our results with the two-sided symmetric case.The reduced model through the one-sided symmetric projection, however, would be the same as obtained from the use of regular form.Results for our proposed one-sided and two-sided regular projection approaches The results show that the additional use of the third transfer function improves the error of both the one-sided projection approach and the two-sided approach.Although the improvement is not significant in the one-sided projection, we gain basically one order of accuracy for the two-sided approach.We believe that this might be due to the choice of interpolation points.Since in the twosided approach, we are matching derivatives as well for the same choice of interpolation points, the effect of interpolating the third transfer function is more obvious in the two-sided case.Now we change the system input to () = cos(2(/10) + 1)/2 and use the same reduced models as before for the = 2 case to obtain the results are shown in Figs.5a and b.Unlike trajectory based methods for model order reduction, the results show that the approximation quality of the projection based reduced models is not effected by the variation in the control input.Next we consider interpolation of the first five subsystems.The size of the reduced system is now = 21 (which should be ≤ 25).The results are shown in Figs.6a and b.Although we are using the same set of interpolation points as those used in Fig. 4, here we get a significant improvement of the one-sided and two-sided approximation errors.This is because we are now interpolating the first five subsystems in the regular form.
The results clearly show that the additional interpolation of the third subsystem improves the approximation error significantly in both one-sided and two-sided projection approaches.The quality of the reduced model varies with the choice of the interpolation points.We used IRKA on the linear part of the system to select and fix a choice of interpolation points.

Conclusion
We extended the orthogonal and oblique projection framework for model reduction of quadratic-bilinear systems to multimoment matching of higher subsystems.For this, we derive the regular form of the multivariate transfer functions associated with the quadratic-bilinear system.The structure of the regular multivariate transfer functions is simpler as compared to the symmetric form of the transfer functions.Multimoment matching of the regular form is therefore easy to ensure for higher subsystems.The choice of interpolation points is an important issue.We selected the interpolation points so that the basis vectors can be reused for other basis vectors.An important future work would be to improve the choice of interpolation points.