## 3 Deriving Input-Output Mappings for the Two Classes of Systems  
### The Variational Equation Approach  
Consider a nonlinear system \(\mathbf{\Sigma}_{\rm N}\) as defined in (1). Assume that the nonlinear functions \(\mathbf{f},\mathbf{g}:\mathbb{R}^{n}\to\mathbb{R}^{n}\) are analytic in \(\mathbf{x}\). Without loss of generality we also enforce that: \(\mathbf{x}(0)=\mathbf{0},\ \mathbf{f}(\mathbf{x}(0))=\mathbf{f}(0)=\mathbf{0}\). Next, compute an approximate bilinear model \(\mathbf{\Sigma}_{\rm B}\) via Carleman linearization presented in Section 2.1.1 and an exact quadratic bilinear model \(\mathbf{\Sigma}_{\rm QB}\) via the reformulation procedure presented in Section 2.2.1. The idea is to decouple these two classes of systems and write them equivalently as an infinite linear time-varying representation. This approach (see [37; 38]) studies the response of the systems to control inputs of the form \(\mathbf{u}(t)=\alpha\hat{\mathbf{u}}(t)\) (\(\alpha\) is an arbitrary scalar, and \(\hat{\mathbf{u}}(t)\) is a unitary signal). Since the dynamical systems under considerations contain analytic non-linearities, it follows that the solution \(\mathbf{x}(t)\) also exhibits an analytic representation and thus can be written as an infinite power series expansion in the parameter \(\alpha\), e.g.,  
\[\mathbf{x}(t)=\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t)=\alpha\mathbf{x} _{1}(t)+\alpha^{2}\mathbf{x}_{2}(t)+\alpha^{2}\mathbf{x}_{3}(t)+\ldots \tag{19}\]  
We proceed with the bilinear model obtained via Carleman linearizaton, i.e., as introduced in (2)). By substituting (19) and the control input \(\mathbf{u}(t)=\alpha\hat{\mathbf{u}}(t)\) onto the differential equation in (2), it follows that:  
\[\mathbf{E}\sum_{k=1}^{\infty}\alpha^{k}\hat{\mathbf{x}}_{k}(t) =\mathbf{A}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t)+\mathbf{ N}\Big{(}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{x}_{k}(t)\Big{)}\alpha\hat{ \mathbf{u}}(t)+\alpha\mathbf{B}\hat{\mathbf{u}}(t)\] \[=\Big{(}\sum_{k=1}^{\infty}\alpha^{k}\mathbf{A}\mathbf{x}_{k}(t) \Big{)}+\Big{(}\sum_{k=2}^{\infty}\alpha^{k}\mathbf{N}\mathbf{x}_{k-1}(t) \mathbf{u}(t)\Big{)}+\alpha\mathbf{B}\hat{\mathbf{u}}(t).\]  
Now, since this differential equation and the equation for the initial state must hold for any value of \(\alpha\), coefficients of same powers of \(\alpha\) can be equated. By selecting the terms corresponding to \(\alpha^{k}\) from both sides, we start with the equality \(\mathbf{E}\mathbf{x}_{1}(t)=\mathbf{A}\mathbf{x}_{1}(t)+\mathbf{B}\hat{ \mathbf{u}}(t)\) (for \(k=1\)). In general, by collecting all terms corresponding to \(\alpha^{k}\) from both sides, for \(k>1\), we obtain equalities, such as \(\mathbf{E}\mathbf{x}_{k}(t)=\mathbf{A}\mathbf{x}_{k}(t)+\mathbf{N}\mathbf{x}_{ k-1}(t)\hat{\mathbf{u}}(t)\). Each equality can be rewritten as \(\mathbf{E}\mathbf{x}_{k}(t)=\mathbf{A}\mathbf{x}_{k}(t)+\mathbf{B}_{k}(t)\hat {\mathbf{u}}(t)\), by introducing \(\mathbf{B}_{k}(t):=\mathbf{N}\mathbf{x}_{k-1}(t)\mathbf{u}(t),\ k>1\). Hence, a linear differential equation with time-varying coefficients (the vector \(\mathbf{B}_{k}(t)\)). The next step is to decouple the differential equation corresponding to the BTI system into infinite many time-varying linear differential equations. We can, hence, rewrite the initial dynamics equivalently as follows:  
\[\mathbf{E}\mathbf{x}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{N}\mathbf{x}(t) \mathbf{u}(t)\mathbf{B}\hat{\mathbf{u}}(t) \Leftrightarrow \begin{cases}\mathbf{E}\mathbf{x}_{1}(t)&=\mathbf{A}\mathbf{x}_{1}(t )+\mathbf{B}\hat{\mathbf{u}}(t),\\ \mathbf{E}\mathbf{x}_{2}(t)&=\mathbf{A}\mathbf{x}_{2}(t)+\mathbf{N}\mathbf{x} _{1}(t)\hat{\mathbf{u}}(t),\\ &\vdots\\ \mathbf{E}\mathbf{x}_{k}(t)&=\mathbf{A}\mathbf{x}_{k}(t)+\mathbf{N}\mathbf{x} _{k-1}(t)\hat{\mathbf{u}}(t),\\ &\vdots\end{cases} \]  
Denote with \(\mathbf{y}_{k}(t)=\mathbf{C}\mathbf{x}_{k}(t),\ \forall\ k\geqslant 1\). The initial BTI system is written equivalently as infinitely many coupled linear subsystems.