## 2 Transforming the Original Dynamics to Quadratic or to Bilinear Format
### Bilinear Systems
#### 2.1.1 Carleman's Linearization  
In what follows, we introduce a well-established numerical tool that enforces approximating a nonlinear system with analytical nonlinearities be means of a bilinear system. The technique is commonly known as Carleman's linearization and was originally introduced in [32]. Since this is an approximation technique, it follows that the bilinear system obtained is clearly not equivalent to the initial one (in terms of the input-output response). One main property of this method is that, as the dimension of the bilinear system increases, the original nonlinear system is better approximated. Of course, in many cases, one cannot overly increase the dimensionality, since this results in challenges, such as low computational time or depleted memory storage. The Kronecker product of the internal state variable \(\mathbf{x}=[x_{1}\ \cdots\ x_{n}]^{T}\) with itself is defined as follows:  
\[\mathbf{x}\otimes\mathbf{x}=[x_{1}^{2}\ x_{1}x_{2}\ \ \ldots\ \ x_{1}x_{n}\ \ \ldots\ \ x_{n}^{2}]^{T}\in\mathbb{R}^{n^{2}}. \tag{3}\]  
Denote the composition of \(k\) Kronecker products of the same vector \(\mathbf{x}\) with itself  
Now, by using a truncated version of the Maclaurin series of the functions by keeping only the first N terms for both series, where \(N\) represents the truncation index, we can write explicitly:  
\[ \begin{cases}\mathbf{f}(\mathbf{x})=\sum_{k=1}^{N}\mathbf{F}_{k}\mathbf{x}^{(N) }=\mathbf{F}_{1}\mathbf{x}+\mathbf{F}_{2}\mathbf{x}^{(2)}+\ldots+\mathbf{F}_{N }\mathbf{x}^{(N)},\\ \mathbf{g}(\mathbf{x})=\sum_{k=0}^{N-1}\mathbf{G}_{k}\mathbf{x}^{(k)}= \mathbf{G}_{0}+\mathbf{G}_{1}\mathbf{x}+\ldots+\mathbf{G}_{N-1}\mathbf{x}^{(N -1)}.\end{cases}  \tag{6}\]  
In the first equality in (6), we can neglect the term \(\mathbf{F}_{0}\), without the loss of generality. This is usually chosen to be \(\mathbf{0}\) in the case of non-existing forcing terms (and if indeed it is non-zero, it can be incorporated in the state vector by means of shifting). By substituting (6) into the state equation in (1), we write  
\[\dot{\mathbf{x}}(t)=\sum_{k=1}^{N}\mathbf{F}_{k}\mathbf{x}^{(k)}(t)+\sum_{k=0 }^{N-1}\mathbf{G}_{k}\mathbf{x}^{(k)}(t)\mathbf{u}(t). \tag{7}\]  
Hence, it follows that the left hand side of the differential equation in (7) contains the state variable vector \(\mathbf{x}(t)\), while the right hand side contains various powers of it \((\mathbf{x}^{(2)}(t),\ldots,\mathbf{x}^{(N)}(t))\). The strategy is to modify the state variable by adding higher-order terms that were missing; in this way, the dimension of the state vector increases. First, we compute the derivative of \(\mathbf{x}^{(2)}(t)\) explicitly, as:  
\[\dot{\mathbf{x}}^{(2)}(t) =\frac{d\mathbf{x}^{(2)}(t)}{dt}=\frac{d\big{(}\mathbf{x}(t) \otimes\mathbf{x}(t)\big{)}}{dt}=\dot{\mathbf{x}}(t)\otimes\mathbf{x}(t)+ \mathbf{x}(t)\otimes\dot{\mathbf{x}}(t)\] \[=\Big{[}\sum_{k=1}^{N}\mathbf{F}_{k}\mathbf{x}^{(k)}(t)+\sum_{k=0 }^{N-1}\mathbf{G}_{k}\mathbf{x}^{(k)}(t)\mathbf{u}(t)\Big{]}\otimes\mathbf{x}(t)\] \[+\mathbf{x}(t)\otimes\Big{[}\sum_{k=1}^{N}\mathbf{F}_{k}\mathbf{ x}^{(k)}(t)+\sum_{k=0}^{N-1}\mathbf{G}_{k}\mathbf{x}^{(k)}(t)\mathbf{u}(t)\Big{]}\] \[=\sum_{k=1}^{N-1}[\mathbf{F}_{k}\otimes\mathbf{I}_{n}+\mathbf{I} _{n}\otimes\mathbf{F}_{k}]\mathbf{x}^{(k+1)}(t)+\sum_{k=0}^{N-2}[\mathbf{G}_{ k}\otimes\mathbf{I}_{n}+\mathbf{I}_{n}\otimes\mathbf{G}_{k}]\mathbf{x}^{(k+1)}(t) \mathbf{u}(t).\]  
We introduce the following summation of \(j\) terms, each containing \(j-1\) Kronecker products, as follows (for \(j\geqslant 2\) and \(k\geqslant 1\))  
\[\mathbf{F}_{j,k}=\mathbf{F}_{k}\otimes\mathbf{I}_{n}\otimes\ldots\otimes \mathbf{I}_{n}+\mathbf{I}_{n}\otimes\mathbf{F}_{k}\otimes\mathbf{I}_{n} \otimes\ldots\otimes\mathbf{I}_{n}+\ldots+\mathbf{I}_{n}\otimes\mathbf{I}_{n} \otimes\ldots\otimes\mathbf{F}_{k},\]  
where we set \(\mathbf{F}_{1,k}:=\mathbf{F}_{k}\), by definition.