## 2 Sensitivity Analysis of Chaotic Systems Using MSS  
In this work, we consider a non-linear dynamical system governed by the equation  
\[\frac{d\mathbf{u}}{dt}=\mathbf{f}(\mathbf{u},\mathbf{s}), \tag{1}\]  
where \(\mathbf{u}\) denotes the state vector, and \(\mathbf{s}\) a set of design parameters. This is a set of ordinary differential equations (ODEs) that arise after spatial discretisation of partial differential equations (PDEs). The QoI is the time average of some function \(G(\mathbf{u},\mathbf{s})\) of the state variables and design parameters, as in  
\[\overline{G}_{\infty}(\mathbf{s})=\lim_{T\to\infty}\frac{1}{T}\int_{0}^{T}G( \mathbf{u},\mathbf{s})dt. \tag{2}\]  
An example QoI from the field of aerodynamics is the time-average of the lift or the drag coefficient of an airfoil. The conventional adjoint sensitivity analysis to determine \(\frac{d\overline{G}_{\infty}(\mathbf{s})}{d\mathbf{s}}\) requires the formation of the augmented Lagrangian  
\[\Lambda=\frac{1}{T}\int_{0}^{T}G(\mathbf{u},\mathbf{s})+\lambda^{T}\left( \frac{d\mathbf{u}}{dt}-\mathbf{f}(\mathbf{u},\mathbf{s})\right)dt+\mu^{T} \left(\mathbf{u}(0)-\mathbf{u}_{0}\right). \tag{3}\]After differentiation with respect to \(\mathbf{s}\) and integration by parts, we arrive at the adjoint system  
\[\frac{d\lambda}{dt}=-\Big{(}\frac{\partial\mathbf{f}}{\partial\mathbf{u}}\Big{)} ^{T}\lambda+\frac{\partial G(\mathbf{u},\mathbf{s})}{\partial\mathbf{s}}, \tag{4}\]  
which is integrated backwards in time with terminal condition \(\lambda(T)=0\). As mentioned in the Introduction, if the system in Equation (1) is chaotic and has at least one positive Lyapunov exponent, then marching Equation (4) backwards will result in an exponentially growing adjoint state vector \(\lambda\). To resolve the problem of diverging trajectories for \(\mathbf{s}\) and \(\mathbf{s}+\delta\mathbf{s}\), we formulate a minimization problem that keeps the two trajectories close in phase space. To achieve this, we relax the requirement of fixed initial condition for the two nearby values of \(\mathbf{s}\). However, for ergodic systems, this change does not affect the time-average values and their sensitivities. Below, we state the minimization problem for a single parameter \(s\), but an adjoint version for multiple parameters can be also formulated. \[\min_{\mathbf{v}(t_{i}^{+})}\sum_{i=0}^{K}\|\mathbf{v}(t_{i}^{+}) \|_{2}^{2}, \tag{5a}\] \[\text{subject to}\ \mathbf{v}(t_{i}^{+})=\mathbf{v}(t_{i}^{-}),\] (5b) \[\frac{d\mathbf{v}}{dt}-\frac{\partial\mathbf{f}}{\partial \mathbf{u}}\mathbf{v}-\frac{\partial\mathbf{f}}{\partial\mathbf{s}}-\eta \mathbf{f}=0,\] (5c) \[\left(\mathbf{f}(\mathbf{u}(t),s),\mathbf{v}(t)\right)=0,\] (5d) \[\text{where}\ t_{i}<t<t_{i+1}\ \left(i=0,...,K-1\right), \tag{5e}\]  
where  
\[ \begin{split}\mathbf{v}(t)&=\lim_{\delta\to 0} \frac{\mathbf{u}(\tau(t);s+ds)-\mathbf{u}(t;s)}{\delta s}\\ \eta(t)&=\frac{d}{ds}\left(\frac{d\tau(t)}{dt}-1 \right).\end{split}  \tag{6}\]  
The objective to be minimized, given in Equation (5a), is the \(L_{2}\) norm of the distance between the two trajectories in phase space. This norm is computed using the distance at \(K+1\) discrete time points \(t_{0},\ldots,t_{K}\). Equation (5b) enforces the continuity of \(\mathbf{v}\) across consecutive segments, Equation (5c) is the linearized version of Equation (1) around \(s\), and Equation (5d) enforces the direction \(\mathbf{v}\) to be normal to \(\mathbf{f}\). Finally, \(\eta\) is a time-dilation term between the reference and the shadowing trajectory. We briefly describe the solution of the minimization problem in Equation (5), and more details can be found in Reference [21; 22]. We reformulate Equation (5) as  
\[\text{Minimize}\ \ \frac{1}{2}\sum_{i=0}^{K}\|\mathbf{v}_{i}\|_{2}^{2}, \tag{7a}\] \[\text{subject to}\ \ A\underline{\mathbf{v}}=\underline{\mathbf{b}}, \tag{7b}\]  
where matrix \(A\) is  
\[A= \begin{bmatrix}-\Phi_{1}&I&&\\ &-\Phi_{2}&I&&\\ &&\ddots&\ddots&\\ &&&-\Phi_{K}&I\end{bmatrix} \ \ \ \underline{\mathbf{v}}=