Multiple Shooting Shadowing for Sensitivity Analysis of Chaotic Dynamical Systems

TL;DR

This paper introduces multiple shooting shadowing (MSS), a more efficient method for sensitivity analysis in chaotic systems, reducing computational costs compared to the least squares shadowing (LSS) approach.

Contribution

The paper proposes MSS, an improved shadowing method that enhances computational efficiency and reduces memory usage for sensitivity analysis in chaotic dynamical systems.

Findings

MSS converges faster than LSS.

MSS requires less memory and computational time.

Applicable to chaotic vortex shedding and turbulence simulations.

Abstract

Sensitivity analysis methods are important tools for research and design with simulations. Many important simulations exhibit chaotic dynamics, including scale-resolving turbulent fluid flow simulations. Unfortunately, conventional sensitivity analysis methods are unable to compute useful gradient information for long-time-averaged quantities in chaotic dynamical systems. Sensitivity analysis with least squares shadowing (LSS) can compute useful gradient information for a number of chaotic systems, including simulations of chaotic vortex shedding and homogeneous isotropic turbulence. However, this gradient information comes at a very high computational cost. This paper presents multiple shooting shadowing (MSS), a more computationally efficient shadowing approach than the original LSS approach. Through an analysis of the convergence rate of MSS, it is shown that MSS can have lower memory usage and run time than LSS.

Figures

Tables

Table 1. Table 1: Variable definitions for Dowell’s buckled plate [ 52 ] . Note that D = E ​ h 3 / 12 ​ ( 1 − ν 2 ) 𝐷 𝐸 superscript ℎ 3 12 1 superscript 𝜈 2 D=Eh^{3}/12(1-\nu^{2}) and N x = ( E ​ h / 2 ​ a ) ​ ∫ 0 a ( ∂ w / ∂ x ) 2 ​ 𝑑 x subscript 𝑁 𝑥 𝐸 ℎ 2 𝑎 superscript subscript 0 𝑎 superscript 𝑤 𝑥 2 differential-d 𝑥 N_{x}=(Eh/2a)\int_{0}^{a}(\partial w/\partial x)^{2}\ dx .

Variable	Definition
$w$	plate transverse deflection
$x$	streamwise spatial coordinate
$t$	time
$a$	plate length
$h$	plate thickness
$m$	mass per unit length of the plate
$N_x^E$	externally applied in-plane load (positive in tension)
$E$	modulus of elasticity
$\nu$	Poisson ratio of the plate material
$D$	plate bending stiffness
$N_x$	the tension created by the stretching of the plate due to bending
${\rho }_{\infty }$	fluid mass density
$U_{\infty }$	flow velocity
$M$	flow mach number
$\Delta p$	static pressure difference across the plate

Equations

$$\frac{du}{dt}=f(u;s),$$

$$J(t;s)=J(u(t;s))$$

$$\bar{J}(s)=\frac{1}{T_1-T_0}{\int }_{T_0}^{T_1}J(t;s)dt$$

$$\frac{du}{dt}$$

$$\frac{d{u}^{\prime }}{dt}$$

$$\frac{d\bar{J}}{ds}\approx \frac{1}{\Delta T}{\int }_{T_0}^{T_0+\Delta T}\frac{J(u^{\prime })-J(u)}{\delta s}dt$$

$$u^{\prime }(t)-u(t)\sim e^{{\Lambda }^{max}t}$$

$$\frac{d\bar{J}}{ds}=\underset{\delta s\to 0}{\lim }\underset{\Delta T\to \infty }{\lim }\frac{1}{\Delta T}{\int }_{T_0}^{T_0+\Delta T}\frac{J(u^{\prime })-J(u)}{\delta s}dt\ne \underset{\Delta T\to \infty }{\lim }\frac{1}{\Delta T}{\int }_{T_0}^{T_0+\Delta T}\underset{\delta s\to 0}{\lim }\frac{J(u^{\prime })-J(u)}{\delta s}dt$$

$$\frac{d\bar{J}}{ds}=\underset{\delta s\to 0}{\lim }\underset{\Delta T\to \infty }{\lim }\frac{1}{\Delta T}{\int }_{T_0}^{T_0+\Delta T}\frac{J(u^{\prime })-J(u)}{\delta s}dt\ne \underset{\Delta T\to \infty }{\lim }\frac{1}{\Delta T}{\int }_{T_0}^{T_0+\Delta T}\left(\frac{\partial J}{\partial u}v(t)+\frac{\partial J}{\partial s}\right)dt$$

$$\underset{u,\tau }{\min }\frac{1}{2}{\int }_{T_0}^{T_1}\Vert u^{\prime }(\tau (t)-u(t){\Vert }^2+\alpha \Vert 1-d\tau /dt{\Vert }^{2}dts.t.\frac{d{u}^{\prime }}{d\tau }=f(u^{\prime };s+\delta s),T_0\le t\le T_1$$

$$\underset{v,\eta }{\min }\frac{1}{2}{\int }_{T_0}^{T_1}\Vert v(t){\Vert }^2+{\alpha }^2\Vert \eta {\Vert }^{2}dt,\text{s.t.}\frac{dv}{dt}=\frac{\partial f}{\partial u}v+\frac{\partial f}{\partial s}+\eta f,T_0\le t\le T_1$$

$$\frac{\partial w}{\partial t}$$

$${\alpha }^2\eta =⟨f,w⟩$$

$$\frac{dv}{dt}$$

$$\underset{v(t_i)}{\min }\sum _{i=0}^K\Vert v(t_i^{+}){\Vert }_2^2$$

$$\text{s.t.}v(t_i^{+})=v(t_i^{-}),i=1,2,...,K$$

$$\frac{dv}{dt}=\frac{\partial f}{\partial u}v+\frac{\partial f}{\partial s}+\eta f,t_0<t<t_K$$

$$⟨f(u(t);s),v(t)⟩=0$$

$$\phi^{t,t}=I\mbox{ for all }t\;,\quad\phi^{t^{\prime},\tau}\cdot\phi^{t,t^{\prime}}=\phi^{t,\tau}\;,\quad\frac{d}{d\tau}\phi^{t,\tau}=\frac{\partial f}{\partial u}\bigg{|}_{\tau}\phi^{t,\tau}\;,\quad\frac{d}{dt}\phi^{*\,t,\tau}=-\frac{\partial f}{\partial u}\bigg{|}_{t}^{*}\phi^{*\,t,\tau}\;,$$

$$v(t)=\phi^{t_{i},t}v(t_{i})+\left(\int_{t_{i}}^{t}\eta_{i}(\tau)\,d\tau\right)f(u(t);s)+\int_{t_{i}}^{t}\phi^{\tau,t}\frac{\partial f}{\partial s}\bigg{|}_{\tau}\,d\tau\;,\quad t_{i}\leq t<t_{i+1}\;,$$

$${\int }_{t_i}^t{\eta }_i(\tau )d\tau =-\frac{⟨v^{\prime }(t),f(u(t);s)⟩}{\Vert f(u(t);s){\Vert }_2^2}$$

$$v^{\prime}(t)=\phi^{t_{i},t}\textbf{v}_{i}+\int_{t_{i-1}}^{t_{i}}\phi^{\tau,t}\frac{\partial f}{\partial s}\bigg{|}_{\tau}\,d\tau$$

$$v(t)=v^{\prime }(t)-\frac{⟨v^{\prime }(t),f(u(t);s)⟩}{\Vert f(u(t);s){\Vert }_2^2}f(u(t);s)t_i\le t<t_{i+1},$$

$$v(t)=P_t{v}^{\prime }(t)\equiv v^{\prime }(t)-\frac{⟨f(u(t);s),v^{\prime }(t)⟩}{\Vert f(u(t);s){\Vert }_2^2}f(u(t);s)$$

$$\underset{{\text{v}}_i}{\min }\sum _{i=0}^K\Vert {\text{v}}_i{\Vert }_2^2$$

$$\displaystyle\text{s.t.}\quad\textbf{v}_{i+1}=P_{t_{i+1}}\phi^{t_{i},t_{i+1}}{\bf v}_{i}+P_{t_{i+1}}\int_{t_{i}}^{t_{i+1}}\phi^{\tau,t}\frac{\partial f}{\partial s}\bigg{|}_{\tau}\,d\tau$$

$${\text{v}}_{i+1}={\Phi }_{i+1}{\text{v}}_i-{\mathbf{b}}_{i+1}$$

$${\bf b}_{i}=-P_{t_{i}}\int_{t_{i-1}}^{t_{i}}\phi^{\tau,t}\frac{\partial f}{\partial s}\bigg{|}_{\tau}\,d\tau$$

$$\text{uuline}\mathbf{B}\underline{\mathbf{v}}=\underline{\mathbf{b}}$$

$$\uuline{{\bf B}}=\left(\begin{array}[]{ccccc}\Phi_{1}&-I&&&\\ &\Phi_{2}&-I&&\\ &&\ddots&\ddots&\\ &&&\Phi_{K}&-I\end{array}\right),\>\underline{{\bf v}}=\left(\begin{array}[]{c}{\bf v}_{0}\\ {\bf v}_{1}\\ \vdots\\ {\bf v}_{K}\end{array}\right),\>\underline{{\bf b}}=\left(\begin{array}[]{c}{\bf b}_{1}\\ {\bf b}_{2}\\ \vdots\\ {\bf b}_{K}\end{array}\right)$$

$$Q\equiv \sum _{i=0}^K\Vert {\text{v}}_i{\Vert }_2^2={\underline{\mathbf{v}}}^T\underline{\mathbf{v}}$$

$$\Lambda ={\underline{\mathbf{v}}}^T\underline{\mathbf{v}}+{\underline{\mathbf{w}}}^T(\underline{\mathbf{b}}-\text{uuline}\mathbf{B}\underline{\mathbf{v}})$$

$$\frac{\partial \Lambda }{\partial \underline{\mathbf{v}}}=$$

$$\frac{\partial \Lambda }{\partial \underline{\mathbf{w}}}=$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Computational methods for sensitivity analysis are invaluable tools for research and design in many engineering and scientific fields. These methods compute derivatives of outputs with respect to inputs in computer simulations. In applications with large amounts of parameters and only a few important outputs, adjoint-based sensitivity analysis is especially efficient [1]. In aircraft design, for example, the number of geometric parameters that define the outer mold line is typically very large, but engineers may only be interested in a few outputs, such as the lift-to-drag ratio. As a result, the adjoint method of sensitivity analysis has proven to be very successful for aircraft design with gradient-based optimization [2, 3, 4]. The adjoint method has also been an essential tool for adaptive grid methods for solving partial differential equations (PDE’s) [5], error estimation [6], and flow control problems [7]. Finally, some techniques for uncertainty quantification can benefit immensely from sensitivity information [8, 9].

Unfortunately, conventional sensitivity analysis methods, including the adjoint method, can break down when applied to chaotic systems. This occurs for sensitivities of long-time-averaged quantities of interest to design inputs. In this context “long-time” refers to time averaging horizons much larger than the physical time scales associated with the chaotic system being considered. This is problematic, as many key scientific and engineering quantities of interest in chaotic systems are long-time-averaged quantities, such as the time-averaged lift or drag coefficient of a flight vehicle in a high-lift configuration or the average heat transfer to a turbine blade.

To carry out efficient design and analysis of chaotic systems and fluid flows, a new sensitivity analysis method is needed. One promising new approach is Least Squares Shadowing (LSS) [10, 11, 12, 13]. The most common implementation of LSS in the literature is called transcription LSS. Transcription LSS involves solving a globally coupled space-time problem, which can be computationally intensive [10]. The study of LSS for chaotic vortex shedding by Blonigan et al. [14] shows that transcription LSS is very costly in memory usage and operation count for a relatively small simulation. Similar issues with computational cost were encountered in a study of transcription LSS for a direct numerical simulation (DNS) of homogeneous isotropic turbulence [15].

This paper presents an alternative way to pose the minimization statement to compute the shadowing direction or the adjoint shadowing direction. This formulation, multiple shooting shadowing (MSS), addresses the high computational cost of LSS. It is shown that MSS can reduce the memory requirements and the run time required to compute sensitivities of chaotic systems, making LSS sensitivity analysis more tractable for large chaotic dynamical systems such as turbulent fluid flow simulations. The convergence properties of MSS are presented in great detail, along with some approaches to control the convergence rate of MSS.

This paper is organized as follows: section 2 discusses how conventional sensitivity analysis approaches break down for chaotic systems and the past work done to avoid this break down. Section 3 provides an overview of the formulation of transcription LSS. Section 4 introduces the MSS minimization statement and section 5 discusses how MSS can be implemented. Next, section 6 shows MSS results for a chaotic dynamical system and a chaotic partial differential equation (PDE). Section 7 discusses the convergence properties of MSS. Finally, section 8 summarizes this thesis and discusses some future research directions.

2 Sensitivity Analysis of a Chaotic dynamical System

The first question to be asked is if sensitivities are in fact well defined for chaotic systems. It is believed that many, but not every, chaotic system has differentiable time averaged quantities $\overline{J}$. Specifically, chaotic systems classified as uniformly hyperbolic or quasi-hyperbolic have differentiable time averaged quantities, but non-hyperbolic systems do not. These classes are discussed in greater detail below.

A uniformly hyperbolic attractor is a strange attractor with a tangent space that can be decomposed into stable, neutrally stable, and unstable subspaces at every point in phase space [16]. In other words, the Lyapunov covariant vectors make up a basis for phase space at all points on an attractor. Ruelle’s linear response theorem states that hyperbolic attractors have mean quantities that respond differentiably to small perturbations to their parameters [17]. Therefore, sensitivities are well defined for chaotic systems with hyperbolic attractors. A well studied example of a hyperbolic attractor is the Plykin attractor [18]. The equations governing the Plykin attractor were designed to have hyperbolic properties, which are rare in practice.

Although uniformly hyperbolic attractors are rare, many important properties of hyperbolic systems, including Ruelle’s linear response theorem, can also be shown to hold for the far more common non-uniformly hyperbolic or quasi-hyperbolic attractors [17, 19, 20]. One example of a quasi-hyperbolic attractor is the Lorenz attractor [21]. At the origin of phase space, the Lyapunov covariant vectors for the positive and negative exponent are parallel, so hyperbolicity does not apply. However, this point is an unstable saddle point and almost all phase space trajectories do not pass through it. Because of this the Lorenz attractor appears to have the properties of a hyperbolic attractor, most importantly differentiable mean quantities [22, 10].

Other chaotic dynamical systems have non-hyperbolic attractors. In these non-hyperbolic systems the time averaged quantities are usually not differentiable or even continuous as the parameters vary. In fact, long-time-averages for non-hyperbolic systems may have nontrivial dependence on the initial condition (i.e. the system is not ergodic), which leads to time averaged quantities that are not well-defined.

Fortunately, Gallavotti and Cohen’s chaotic hypothesis conjectures that larger systems behave more like hyperbolic systems than non-hyperbolic systems [23, 24]. That is, larger systems should have differentiable infinite time-averaged quantities. Additionally, a study by Albers and Sprott found that larger chaotic systems tend to have smoothly varying topology changes in the attractor as system parameters are varying [25]. Long-time-averaged quantities do not necessarily vary smoothly across sudden topology changes like bifurcations. Therefore, the chaotic hypothesis and the work by Albers and Sprott suggest there are well defined sensitivities to be computed for a large range of chaotic systems, especially if these systems have a large number of degrees of freedom (DoF). This is encouraging considering the large numbers of DoF’s in simulations such as those of chaotic and turbulent fluid flows.

Additionally, there is some evidence that the chaotic hypothesis applies to simulations of turbulent fluid flows. Grid convergence studies have been done in many cases to ensure that the discretization of the governing equations is sufficiently detailed. For example, Kim et al. used a direct numerical simulation (DNS) to compute turbulent statistics such as the mean velocity profile of a turbulent channel flow with a coarse and a fine spatial discretization [26] to check if their fine discretization was fine enough. The statistics were the same for the coarse and fine discretizations, which shows that long-time-averaged quantities of the DNS respond smoothly to perturbations in the spatial discretization, as predicted by the chaotic hypothesis. Similar results of grid convergence studies for other DNS and Large Eddy Simulation (LES) results [27, 28] also support the chaotic hypothesis, making it very likely that sensitivities of long-time-averaged quantities are well defined for high-fidelity turbulent flow simulations.

2.1 Conventional Sensitivity Analysis

Many time dependent simulations can be interpreted as dynamical systems:

[TABLE]

where $u$ is a length $n$ vector representing the system state and $s$ is a set of design parameters. For a computational fluid dynamics (CFD) simulation, equation (1) represents the discretized Navier-Stokes equations and the state variable $u$ is a vector containing the conserved quantities. The design parameters $s$ could include geometric parameters for a wing or flow parameters such as the Reynolds number.

When conducting design studies, engineers are typically interested in minimizing some objective function, $J$:

[TABLE]

One example of $J$ is the instantaneous drag on an airfoil. For unsteady simulations, time-averaged objective functions are often of interest:

[TABLE]

Sensitivity analysis seeks to compute the derivative of $\bar{J}$ with respect to the design parameters $s$. Traditionally, sensitivity analysis is conducted by solving the following initial value problems

[TABLE]

where $u(t)$ is a reference solution and $u^{\prime}(t)$ is solution corresponding to a perturbation $\delta s$ to some design parameter. Sensitivities can be estimated as:

[TABLE]

However, this initial value problem is very poorly conditioned when the system is chaotic [11]. The “butterfly effect” ensures that $u(t)$ and $u^{\prime}(t)$ will become decorrelated after some time. This is because when $u^{\prime}(t)-u(t)$ is infinitesimal:

[TABLE]

where $\Lambda^{max}$ is the largest Lyapunov exponent of the system [29] (see A for one way to compute this). Since chaotic dynamical systems have at least one positive Lyapunov exponent, $u^{\prime}(t)-u(t)$ grows exponentially. Therefore, small perturbations to the initial condition or parameter $s$ will grow relatively large after a time $\mathcal{O}(1/\Lambda^{max})$.

Equation (3) also implies that the linearization $v=\lim_{\delta s\to 0}\frac{u^{\prime}(t)-u(t)}{\delta s}$ will grow exponentially. This growth of $v(t)$ is the cause of the following inequality

[TABLE]

or

[TABLE]

Therefore, sensitivity analysis formulated as an initial value problem does not compute useful sensitivities for chaotic dynamical systems. These same issues arise for conventional adjoint sensitivity analysis as well, as shown by Lea et al. for the Lorenz 63 equation [22].

2.2 Chaotic Sensitivity Analysis

Prior work in sensitivity analysis of long-time-averages in chaotic dynamical systems and fluid flows has been done mostly by the climatological and meteorological community. This work includes the ensemble-adjoint method proposed by Lea et al. [22]. This method has been applied to the Lorenz 63 system and an ocean circulation model [30]. Eyink et al. then went on to generalize the ensemble-adjoint method [19]. The ensemble-adjoint method involves averaging over a large number of adjoint calculations for different segments of a time horizon (or time horizons). It was found by Eyink et al. that the sample mean of sensitivities computed with the ensemble adjoint approach for the Lorenz 63 system converges slower than $N^{-0.5}$, where $N$ is the number of samples, making it less computationally efficient than a naive Monte-Carlo approach [19]. Recently, Ashley and Hicken explored using the ensemble adjoint approach for gradient-based optimization of a chaotic system, but encountered similar issues and observed similarly slow convergence as Eyink et al. [31].

Another idea is the Fokker-Planck adjoint approach for climate sensitivity analysis [32]. This approach involves finding an invariant measure or stationary density which satisfies a Fokker-Planck equation to model the long-time-averaged dynamics, called the “climate” of the system in much of the literature. The adjoint of this Fokker-Planck equation is then used to compute derivatives with respect to long-time-averaged quantities. Fokker-Planck methods typically require discretizing phase space, which limits the method to fairly low dimensional systems. Also many variants of the method require adding diffusion into the system, potentially making the computed sensitivities inaccurate. A Fokker-Planck method proposed by the author and Wang eliminates the need for much of this additional diffusion, but requires a discretization of a manifold approximating the strange attractor, which poses challenges for higher-dimensional, less well understood strange attractors [33].

An analysis based on the Fokker-Planck equation produces the Fluctuation Dissipation Theorem (FDT) [34, 35]. For conservative and nearly conservative dynamical systems, the FDT can be used to accurately compute climate sensitivities [36]. Several improved algorithms based on FDT have since been developed for computing climate sensitivity of non-conservative systems [37, 38, 39]. However, for strongly dissipative systems whose SRB measure [40] deviates strongly from Gaussian, FDT based methods can be inaccurate. Additionally, some of the proposed approaches require computing positive Lyapunov exponents and their corresponding covariant vectors, the current algorithm for which is prohibitively expensive for large systems [38, 41].

3 Transcription Least Squares Shadowing

The poor conditioning of the initial value problem in section 2.1 arises from the problem formulation: $u(t)$ and $u^{\prime}(t)$ will only be close in phase space at $t=T_{0}$ where they share the initial condition $u_{0}$; nothing in the problem formulation requires $u(t)$ and $u^{\prime}(t)$ to be correlated. LSS overcomes the issues of the initial value problem by minimizing the distance in phase space between $u(t)$ and $u^{\prime}(t)$ in a least squares sense [10]. This is done by assuming ergodicity and replacing the initial condition with a regularization, as in Doedel and Friedman’s continuation method for computing heteroclinic orbits [42]. This regularization forces $u^{\prime}(t)$ and $u(t)$ to be as close to one another in phase space as possible for $T_{0}\leq t\leq T_{1}$:

[TABLE]

where $\tau(t)=(1+\delta s\eta t)$ is a time transformation whose purpose is explained in other LSS literature, and $\eta$ called the time dilation term [10].

LSS uses the assumption of ergodicity to convert an initial value problem to a boundary value problem to improve the conditioning of solving for $u^{\prime}(t)$. The solution to this boundary value problem, $u^{\prime}(t)$, called the “shadow trajectory”, has its existence guaranteed by the shadowing lemma for uniformly hyperbolic systems [43, 44, 10].

To compute sensitivities efficiently, equation (5) is linearized [10]:

[TABLE]

where $v(t)=\lim_{\delta s\to 0}\frac{u^{\prime}-u}{\delta s}$ is called the “shadowing direction”. Equation (6) is a linearly constrained least-squares problem with the following Karush-Kuhn-Tucker (KKT) conditions derived using calculus of variations:

[TABLE]

To compute many sensitivities efficiently, one can derive the adjoint equations for equation (6). Adjoint LSS is derived by Wang et al. [10].

The main approach used to numerically solve for the shadowing direction (or the adjoint shadowing direction) is called transcription LSS. Transcription LSS is when equations (7), (8), and (9) are discretized and solved simultaneously for all time $T_{0}<t<T_{1}$ [10]. The KKT system formed by equations (7), (8), and (9) can be reduced to a $mn\times mn$ block tridiagonal linear matrix, where $m$ is the number of time steps and $n$ is the number of degrees of freedom of the system. The KKT system is solved directly for smaller systems and solved using a multigrid in space and time scheme for larger systems [13, 15].

Transcription LSS has been used to compute accurate sensitivities for a computational fluid dynamics (CFD) simulation of chaotic vortex shedding [14] and a DNS of homogeneous isotropic turbulence [15]. Unfortunately, these studies also show that transcription LSS is very costly in memory usage and operation count. Much of this cost is due to the large size of the KKT system. For example, the chaotic vortex shedding simulation has $n=11,090$ degrees of freedom. For the $m=2,000$ time step window considered in the study, the total size of the KKT Schur complement matrix is 22,180,000 rows and roughly 24 GB of input data is required [14]. The system was solved on 2,000 cores for roughly 17 hours to converge the sensitivity to at least 3 decimal places. This high cost provides motivation for the work presented in this paper on the multiple shooting implementation of LSS, called Multiple Shooting Shadowing (MSS).

4 Multiple Shooting Shadowing Formulation

The cost of LSS can be reduced by reformulating the least squares minimization problem (6). More specifically, the objective function in equation (6) can be modified to reduce the number of constraint equations needed to compute a shadowing direction $v(t)$. Instead of minimizing the integral of $v(t)$ over a time horizon as in transcription LSS (6), MSS seeks to minimize $v(t)$ at discrete checkpoints in time.

[TABLE]

Define time segment $i$ as the time span $t_{i-1}\leq t\leq t_{i}$, as shown in figure 1. The MSS minimization problem penalizes large values of $v(t)$ at the end of each time segment. By doing this MSS seeks to find $v(t)$ that does not grow exponentially over each time segment.

Equation (11) shows that MSS has $Kn$ constraint equations for an $n$-DoF system, compared to $mn$ for transcription LSS, where $m$ is the number of time steps used to discretize the time horizon being studied. Since $K$ can be smaller than $m$, the MSS minimization problem can be smaller than that for transcription LSS. In fact, the number of time segments $K$ will almost always be smaller than then total number of time steps $m$. This is because the values of $K$ and $m$ are limited in different ways. If a time segment is too large, $v(t)$ will grow too large between checkpoints and round-off errors will impact the accuracy of MSS. The growth rate of $v(t)$ is at most the largest Lyapunov exponent $\Lambda^{max}$. This means that $\Lambda^{max}$ should be used to determine the minimum value of $K$ for a fixed time horizon. The number of time steps, $m$ is fixed by the desired temporal accuracy and the stability requirements of the time stepping scheme being used. Typically, the time step size needed for stability and accuracy is less than the time needed for round-off errors in the tangent solution to grow to $\mathcal{O}(1)$, so $K<m$. In fact, $K$ can be much smaller than $m$ when the discretization of the governing equation (1) is stiff or the time accuracy requirements are strict. For example, for the Dowell’s plate model presented later in this chapter, $m/K=500$ for a relatively large choice of $K$.

The constraint equation (11) enforces continuity at each checkpoint, which forces $v(t)$ to satisfy the tangent governing equation (12). Equation (13) fixes $\eta(t)$ to eliminate any component of $v(t)$ parallel to $f(u(t);s)$. In other words, equations (12) and (13) combine to form a well-defined algebraic-differential equation, in which $v$ contains the differential variables and $\eta$ is the algebraic variable.

Before proceeding to derive an expression for $\eta$, it is useful to define the tangent propagator $\phi^{t,t^{\prime}}$

[TABLE]

where the superscript $*$ indicates an adjoint or transpose operator. The tangent propagator can be used to write the solution to equations (12) as

[TABLE]

The tangent solution $v(t)$ from equation (15) will satisfy (12) for any $\eta(t)$. For $v(t)$ to also satisfy equation (13), $\int_{t_{i}}^{t}\eta_{i}(\tau)\,d\tau$ must satisfy the following closed form expression derived from equations (15) and (13)

[TABLE]

where

[TABLE]

Equation (16) and the above definition of $v^{\prime}(t)$ can then be substituted into equation (15) to eliminate $\eta(t)$

[TABLE]

Therefore, the tangent solution $v(t)$ can be made to satisfy (13) by using the projection operator $P_{t}$

[TABLE]

To derive the MSS KKT system, equations (15) and (17) are substituted into equations (10) and(11):

[TABLE]

where $\textbf{v}_{i}=v(t_{i}^{+})$.

For a dynamical system with $n$ degrees of freedom, (19) can be written as

[TABLE]

where $\Phi_{i+1}\equiv P_{t_{i+1}}\phi^{t_{i},t_{i+1}}$ is an $n\times n$ matrix called the tangent transition matrix and

[TABLE]

is a $n\times 1$ vector.

All $K$ discretized constraints (20) can be written as a system of equations:

[TABLE]

where

[TABLE]

The block matrix $\uuline{{\bf B}}$ is $Kn$ by $K(n+1)$. The vectors $\underline{\bf v}$, $\underline{\bf w}$, and $\underline{\bf b}$ are length $(K+1)n$, $Kn$, and $Kn$, respectively.

Using this notation, equation (10) can be rewritten as

[TABLE]

Next a Lagrangian is formed for equations (23) and (22).

[TABLE]

This Lagrangian is differentiated with respect to $\underline{\bf v}$ and the Lagrange multiplier $\underline{\bf w}$ to obtain the KKT equations

[TABLE]

Together, the KKT equations form a system of equations that can be solved for $\underline{\bf v}$ and $\underline{\bf w}$:

[TABLE]

As for transcription LSS, the Schur complement of equation (24) is solved instead of the larger full KKT system. The KKT Schur complement for tangent MSS is

[TABLE]

or

[TABLE]

In practice, the matrix $\uuline{{\bf B}}\uuline{{\bf B}}^{T}$ is never formed, since each $\Phi_{i}$ is a dense $n\times n$ matrix and is expensive to compute. Instead, a routine for the product $\uline{\bf y}=\uuline{{\bf B}}\uuline{{\bf B}}^{T}\underline{\bf x}$ is implemented and an iterative method is used to solve (25), as shown in section 5.

Although the matrix $\uuline{{\bf B}}\uuline{{\bf B}}^{T}$ is never formed, its structure reveals some key attributes of MSS. From equations (25) and (26), it can be seen that the Schur complement matrix is a $Kn\times Kn$, symmetric positive definite and block tridiagonal matrix. As stated previously, $K$ is typically much smaller than the total number of time steps $m$. Therefore, the $Kn\times Kn$ MSS KKT Schur complement is potentially much smaller than the $mn\times mn$ KKT Schur complement for transcription LSS.

Additionally, the structure of the MSS KKT Schur complement is independent of the time discretization used. This is different from transcription LSS, for which the structure of the KKT matrix system depends on the time discretization used [10]. The second order Crank-Nicolson scheme considered by Wang et al. [10] leads to a tridiagonal KKT Schur complement, but higher order time discretizations will lead to more bands in the KKT Schur complement matrix. Therefore MSS has a smaller KKT matrix and possibly also a smaller bandwidth than transcription LSS for a given time horizon and time discretization. This can lead to gains in computational efficiency and a potentially simpler implementation, as shown in section 5.

It should be noted that as $K\to\infty$, the MSS minimization problem with uniform time segments becomes equivalent to the transcription LSS minimization problem (6) for $\alpha=0$, where $\alpha$ is the weighting parameter in (6). This is explained in more detail in appendix B.

It should be also be noted that the block structure of equation (25) is similar to the LSS KKT system derived for maps [44]. This is because the operator $\Phi_{i}$ can be thought of as a Jacobian of a nonlinear mapping of $u$ from $t_{i}$ to $t_{i+1}$. For MSS, the nonlinear mapping is solving for $u(t_{i+1})$ from $u(t_{i})$.

As for transcription LSS, the sensitivity of a time-averaged objective function $\bar{J}$ is computed using the following expression [10]

[TABLE]

for an averaging window from $t_{0}$ to $t_{K}$ and $T\equiv t_{K}-t_{0}$. The sensitivity of $\bar{J}$ to $s$ can be written as a function of $v^{\prime}(t)$ from equation (17) (see appendix C for the derivation of the second term)

[TABLE]

where $f_{i+1}=f(u(t_{i+1});s)$ and $J_{i+1}=J(u(t_{i+1});s)$.

4.1 The Filtering Parameter

The MSS minimization problem expressed by equation (10) and (11) can be modified as follows for a better conditioned KKT system

[TABLE]

where $\textbf{v}_{i}=v(t_{i}^{+})$ and $\textbf{w}_{i}=w(t_{i}^{-})$, as shown in figure 1, and $\epsilon$ is called the filtering parameter. Equation (29) shows that non-zero values of $\epsilon$, create a discontinuous jump of $\epsilon\textbf{w}_{i}$ in $v(t)$ at the checkpoints $t_{i}$. In other words, non-zero values of $\epsilon$ relax the continuity constraint (11). This has a major impact on the computational efficiency of MSS that will discussed in more detail in section 7.

To derive the KKT equation for MSS with a non-zero $\epsilon$, start by rewriting equations (28) in the notation of equation (22)

[TABLE]

and (29) as

[TABLE]

where $\underline{\bf w}$ is a length-$Kn$ vector

[TABLE]

The constraint equation (31) can be used to eliminate $\underline{\bf w}$ from equation (30)

[TABLE]

As $Q$ is a quadratic function, it is convex, so the minimum of $Q$ satisfies the first order optimality condition:

[TABLE]

Using equation (31), this expression can be rewritten as

[TABLE]

Together, equations (31) and (33) form a system of equations that can be solved for $\underline{\bf v}$ and $\underline{\bf w}$:

[TABLE]

Finally, the KKT Schur complement for tangent MSS with non-zero $\epsilon$ is

[TABLE]

Note that the left-hand side is identical to that of equation (25) with $\epsilon$ times the identity matrix $\uuline{I}$ added on. In section 7 it is shown that this additional term improves the conditioning of (35), allowing iterative solvers to converge to a solution in less time as $\epsilon$ is increased.

4.2 Adjoint Formulation

When the sensitivities of one time-averaged objective function $\bar{J}$ to many design parameters $s$ are desired, adjoint MSS can be used. Adjoint MSS is the discretely consistent adjoint of tangent MSS and can be expressed as the solution of the following linear system, (see appendix D):

[TABLE]

where $\underline{\bf g}^{T}=({\bf g}_{1}^{T},{\bf g}_{2}^{T},...{\bf g}_{K}^{T},0)$ is a $1\times(K+1)n$ vector and

[TABLE]

The adjoint equations (36) and (37) are derived algorithmically from equations (35) and (27) follwing the discrete adjoint approach as shown in appendix D. Note that the matrix on the left hand side of equation (36) is the same as the one for the tangent MSS KKT Schur complement (35). This is because the matrix $\uuline{{\bf A}}$ is symmetric, and a symmetric matrix is self adjoint.

The sensitivities of $\bar{J}$ can be computed from the adjoint solution $\hat{w}(t)$ with the following expression

[TABLE]

where

[TABLE]

or, in differential equation form

[TABLE]

This can derived by using equation (21), as shown in appendix D.

Interestingly, the adjoint MSS KKT Schur complement (36) is a solution of the following minimization problem

[TABLE]

This is called a ridge regression, or Tikhonov regularization of

[TABLE]

or, equivalently

[TABLE]

this can also be written as

[TABLE]

This shows that equation (42) can be interpreted as solving the conventional adjoint equation (39) in the time horizon $t_{0}<t<t_{K}$, with the initial and terminal conditions

[TABLE]

As $\uuline{{\bf B}}^{T}$ has more rows than columns, equation (42) is over-constrained and can be solved with some kind of regularization. The Tikhonov regularization is easy to analyze because unlike other popular regularizations it has an analytical solution for equation (36), the adjoint MSS KKT Schur complement.

5 Implementation

5.1 Adjoint MSS Algorithm

Like transcription LSS implementations, the MSS implementation is comprised of a solver for the KKT Schur complement and a routine to compute gradients from the tangent or adjoint shadowing direction. For MSS, the KKT Schur complement matrix in equation (35) or (36) is not formed explicitly. Instead, a multiple shooting algorithm is derived for equations (35) or (36). This section presents and discusses the multiple shooting algorithm for carrying out adjoint MSS. This algorithm was used to produce the results presented in this paper. The corresponding tangent MSS algorithm can be found in appendix E.

Adjoint MSS Solver

Inputs: Initial condition for the governing equations $u_{0}$, Spin-up time $t_{0}$, Specified time horizon and checkpoints $t_{0},t_{1},...,t_{K}$, Initial guess for the $Kn\times 1$ vector $\underline{\hat{\textbf{w}}}$, which contains the adjoint variables at checkpoints $1$ to $K$ (default value 0);

Ouputs: Sensitivities $d\bar{J}/ds$

Calls: MATVEC algorithm that computes $\underline{\bf R}=\uuline{{\bf A}}\underline{\hat{\textbf{w}}}+\beta\uuline{{\bf B}}\underline{\bf g}$ where $\underline{\bf R}$ is a $Kn\times 1$ residual vector.

Time integrate the governing equations (1) to compute $u(t)$ for the specified time horizon. Save the objective function $\bar{J}$, it is needed for the right hand side of the linear system (36). 2. 2.

To form the right hand side of the linear system, $\uuline{{\bf B}}\underline{\bf g}$, use the MATVEC algorithm with $\beta=-1$ and $\underline{\bf w}=0$. 3. 3.

Use some iterative algorithm to solve equation (35). To compute the left hand side $\uuline{{\bf A}}\underline{\hat{\textbf{w}}}$, use MATVEC with $\beta=0$. 4. 4.

Compute the sensitivity $d\bar{J}/ds$ using equation (38).

Next, a serial MATVEC algorithm is presented. Note that $t_{i}^{-}$ and $t_{i}^{+}$ refer to the time at checkpoint $i$ in time segments $i-1$ and $i$, respectively. For example, $\hat{w}(t_{i}^{+})$ is the adjoint solution at checkpoint $i$ in time segment $i$.

Serial Adjoint MATVEC Algorithm

Inputs: $\underline{\hat{\textbf{w}}}$, a $Kn\times 1$ vector of the adjoint variables at checkpoints 1 to $K$; $\beta$, a scalar;

Ouputs: $\underline{\bf R}$, a $Kn\times 1$ residual vector

MATVEC computes $\underline{\bf R}=\uuline{{\bf A}}\underline{\hat{\textbf{w}}}+\beta\uuline{{\bf B}}\underline{\bf g}$

For all time segments, compute $\hat{w}(t_{i-1}^{+})$ by integrating $\frac{d\hat{w}}{dt}=-\left(\frac{\partial{f}}{\partial{u}}\right)^{*}\hat{w}+\beta\frac{1}{T}\frac{\partial{J}}{\partial{u}},\;t\in(t_{i-1},t_{i})$ backwards in time with the terminal condition $\hat{w}(t_{i})=P_{t_{i}}\hat{\textbf{w}}_{i}$. 2. 2.

Save $\hat{\textbf{v}}_{i-1}\equiv\hat{w}(t_{i-1}^{+})-\hat{\textbf{w}}_{i-1}$. If $i=1$, save $\hat{\textbf{v}}_{0}\equiv\hat{w}(t_{0}^{+})$. 3. 3.

For all time segments, compute $\hat{v}^{\prime}(t_{i}^{-})$ by integrating $\frac{d\hat{v}^{\prime}}{dt}=\frac{\partial{f}}{\partial{u}}\hat{v}^{\prime},\;t\in(t_{i-1},t_{i}]$ with the initial condition $\hat{v}^{\prime}(t_{i-1})=\hat{\textbf{v}}_{i-1}$. 4. 4.

For all time segments, compute ${\bf R}_{i}=P_{t_{i}}\,\hat{v}^{\prime}(t_{i}^{-})-\hat{\textbf{v}}_{i}+\epsilon\hat{\textbf{w}}_{i}$. If $i=K$, $\hat{\textbf{v}}_{K}=-\hat{\textbf{w}}_{K}$.

Some details of the above algorithms are discussed below.

5.2 Time Integration for MSS

Steps 1 and 3 of the MATVEC algorithm involve solving an adjoint equation for $\hat{w}(t)$ and a tangent equation for $\hat{v}^{\prime}(t)$, respectively. The algorithm presented in this paper should work with any type of time integration scheme except multi-step schemes such as Adams-Bashforth and Adams-Moulton schemes. These scheme require multiple time steps for an initial condition, which the MATVEC algorithm presented in this paper does not support. Also, it is recommended that the scheme used for $\hat{w}(t)$ is the discrete adjoint of the scheme used for $\hat{v}^{\prime}(t)$. Using discretely adjoint time integration schemes for solving the tangent and adjoint in MSS ensures that the matrix represented by the MATVEC algorithm is symmetric.

Additionally, it is recommended that the time integration scheme used for computing $\hat{v}^{\prime}(t)$ is the same as that used for computing the solution to the governing equations $u(t)$. This makes the numerical solutions of $u(t)$ and $\hat{v}^{\prime}(t)$ discretely consistent. If the numerical solutions of $u(t)$, $\hat{v}^{\prime}(t)$, and $\hat{w}(t)$ are all discretely consistent, then verification procedures similar to those used for transcription LSS [14] can be used for MSS. If $u(t)$ and $\hat{v}^{\prime}(t)$ are not discretely consistent, but $\hat{v}^{\prime}(t)$ and $\hat{w}(t)$ are computed by solvers with consistent and stable discretizations, then MSS should still compute accurate sensitivities, but the implementation would be more difficult to verify.

It should be noted that a pair of time integration schemes for $\hat{v}^{\prime}(t)$ and $\hat{w}(t)$ derived or computed using automatic differentiation will be discretely adjoint and discretely consistent with the scheme used to compute the governing equation solution $u(t)$.

5.3 Iterative Solver for MSS

The 3rd step of the Adjoint MSS algorithm calls for some iterative solver for equation (36). The implementation used for the results presented in this thesis uses MINRES, a Krylov subspace method for solving symmetric, sparse matrices [45, 46]. Another Krylov subspace method, conjugate gradient (CG), was tried on a few chaotic systems, but it was found to converge slightly slower than MINRES. Krylov subspace methods such as GMRES are also used by Sanchez and Net in their multiple shooting continuation algorithm [47].

Krylov subspace methods are not the only viable methods to solve the adjoint MSS KKT Schur complement (36). Other iterative solvers based on Multigrid-in-time could also be used, for example, Multigrid Reduction Methods (MGRIT) [48, 49, 50]. MGRIT uses a multiple shooting approach to solve PDEs, with coarse time step solutions on each time segment used to accelerate the convergence of a fine time step solution. Since MGRIT and MSS both incorporate multiple shooting, a promising area for future work is the study of MGRIT approaches for MSS.

5.4 Time-Parallel MSS

The MATVEC algorithm presented in section 5.1 computes the adjoint and then the tangent on each time segment in sequence, but the computations of the adjoint and tangent on each segment can also be done in parallel. The MATVEC algorithm presented below is parallel-in-time and has each time segment assigned to one processor.

Time-Parallel Adjoint MATVEC Algorithm

Inputs: $\underline{\hat{\textbf{w}}}$, a $Kn\times 1$ vector of the adjoint variables at checkpoints 1 to $K$; $\beta$, a scalar;

Ouputs: $\underline{\bf R}$, a $Kn\times 1$ residual vector

Parallel MATVEC computes $\underline{\bf R}=\uuline{{\bf A}}\underline{\hat{\textbf{w}}}+\beta\uuline{{\bf B}}\underline{\bf g}$ on $K$ processors. The $i$th time segment is assigned to processor $i$. Processor $i$ requires the quantities $u(t)$ and $f(u(t);s)$ (or the ability to compute them) for $t\in(t_{i-1},t_{i})$.

To start, processor $i$ has access to $\hat{\textbf{w}}_{i}$.

On processor $i$, do the following:

Compute $\hat{w}(t_{i-1}^{+})$ by integrating $\frac{d\hat{w}}{dt}=-\left(\frac{\partial{f}}{\partial{u}}\right)^{*}\hat{w}+\beta\frac{1}{T}\frac{\partial{J}}{\partial{u}},\;t\in(t_{i-1},t_{i})$ backwards in time with the terminal condition $\hat{w}(t_{i})=P_{t_{i}}\hat{\textbf{w}}_{i}$. Meanwhile, if $i<K$, send $\hat{\textbf{w}}_{i}$ to processor $i+1$. Also, if $i>1$, receive $\hat{\textbf{w}}_{i-1}$ from processor $i-1$. 2. 2.

Compute and save $\hat{\textbf{v}}_{i-1}\equiv\hat{w}(t_{i-1}^{+})-\hat{\textbf{w}}_{i-1}$. If $i=1$, save $\hat{\textbf{v}}_{0}\equiv\hat{w}(t_{0}^{+})$. 3. 3.

Compute $\hat{v}^{\prime}(t_{i}^{-})$ by integrating $\frac{d\hat{v}^{\prime}}{dt}=\frac{\partial{f}}{\partial{u}}\hat{v}^{\prime},\;t\in(t_{i-1},t_{i}]$ with the initial condition $\hat{v}^{\prime}(t_{i-1})=\hat{\textbf{v}}_{i-1}$. If $i=1$, use the initial condition $\hat{v}^{\prime}(t_{0})=\hat{\textbf{v}}_{0}$. Meanwhile, if $i>1$, send $\hat{\textbf{v}}_{i-1}$ to processor $i-1$. Also, if $i<K$, receive $\hat{\textbf{v}}_{i}$ from processor $i+1$. 4. 4.

Compute ${\bf R}_{i}=P_{t_{i}}\,\hat{v}^{\prime}(t_{i}^{-})-\hat{\textbf{v}}_{i}+\epsilon\hat{\textbf{w}}_{i}$. For $i=K$, $\hat{\textbf{v}}_{K}=-\hat{\textbf{w}}_{K}$.

Due to the block tridiagonal structure of equation (36), time-parallel MATVEC has a very straightforward processor communication pattern: processor $i$ only needs to communicate with processors $i-1$ and $i+1$. This means that time-parallel MATVEC will scale very well: the time required for time-parallel MATVEC should be roughly the same as that required to solve an adjoint and then a tangent equation over one time segment. For large dynamical systems or for processors with slow connections there might be some additional time needed to complete communication of $\hat{\textbf{w}}_{i}$ and $\hat{\textbf{v}}_{i}$. For a fixed time segment size, time-parallel MATVEC is much faster than serial MATVEC for large values of $K$ which requires time for $K$ adjoint solves, followed by $K$ tangent solves. For a fixed time horizon $T_{1}-T_{0}$, the run time of the time-parallel MATVEC will decrease if the number of processors increases with $K$.

Note that scaling of a time-parallel MSS solver also depends on the linear solver used for step 3 of the Adjoint MSS solver algorithm. For instance, Krylov subspace solvers like MINRES, GMRES, and Conjugate Gradient require a reduce-all operation to compute dot products of the solution vector $\underline{\hat{\textbf{w}}}$ and residual $\underline{\bf R}$ [45, 46, 51]. Since reduce-all operations involve communication between all processors and one root processor, they have a negative impact on scaling of the linear solver and therefore MSS.

The scaling of time-parallel MSS also depends on the spectral properties of the KKT system which are discussed in section 7. Assume there are one or more time segments per processor. For a fixed time-segment size and a non-zero value of $\epsilon$, the minimum and maximum eigenvalues of the KKT system will be independent of $K$. In this case parallel MSS should scale well if a Krylov solver is used because the condition number will stay the same, so the number of solver iterations should stay the same as well222This assumes the clustering of eigenvalues on the real line does not vary much with $K$. . For a fixed time horizon size $T_{1}-T_{0}$ and a non-zero value of $\epsilon$, the maximum eigenvalue of the KKT system will decrease with $K$, while the minimum eigenvalue will be at least $\epsilon$. In this case, parallel MSS will scale well, since the time to run the MATVEC will decrease and a decrease in condition number means than the KKT system will likely need fewer solver iterations to converge. However, one must be careful of the impact of $\epsilon$ on the accuracy of sensitivities, which is discussed in section 7.3. If $\epsilon=0$ parallel MSS may not scale well as the minimum eigenvalue will decrease as $K$ increases and the KKT system may require more solver iterations to converge.

Finally, it should be mentioned that additional parallelization is possible for governing equations defined over space, like the Navier-Stokes equation. Space-parallel CFD solvers, such as FUN3D, have been in use for a number of years now. Using a space-parallel tangent and adjoint solver in steps 1 and 3 of time-parallel MATVEC would create a space-time-parallel MATVEC. In this case, “processor $i$” would actually refer to a group of processors corresponding to the tangent and adjoint solvers in time segment $i$. If MSS is to be extended to very large chaotic dynamical systems such as LES, a space-time-parallel MATVEC will be necessary to make the time-to-solution reasonable. Therefore, a space-time-parallel MATVEC should be investigated in future studies of MSS.

5.5 Memory requirements of MSS

The adjoint MSS solver algorithm can be implemented to make memory requirements more reasonable than those discussed for the transcription LSS implementation used in FUN3D [14]. This implementation required storing the governing equation solution $u(t)$ and the linearization matrices, including $\partial f/\partial u$ and $\partial f/\partial s$, at all times in the time horizon studied.

For MSS, major savings in memory usage could be realized by using checkpointing for the adjoint solver in step 3 of the MATVEC algorithm [8]. Checkpointing allows one to time integrate the adjoint $\hat{w}(t)$ backwards in time without saving the governing equation solution $u(t)$ and the linearization matrices. Further memory savings can be achieved by solving the tangent equation for $\hat{v}^{\prime}(t)$ in step 3 simultaneously with the governing equations. Such a scheme would only require storage space for $\hat{v}^{\prime}(t)$, $u(t)$, and the linearization matrices required for one time step. Both of these strategies allow for relatively low memory usage when computing the solution of the adjoint MSS KKT Schur complement equation (36).

Computing sensitivities using equation (38) could also be made memory efficient by computing the integral

[TABLE]

on the fly as $\hat{w}(t)$ is solved, instead of solving for $\hat{w}(t)$, storing it, then computing the sensitivities of $\bar{J}$ to different input parameters $s$.

MSS could be implemented by recomputing $u(t)$ on the fly while solving for $\hat{w}(t)$ and $\hat{v}^{\prime}(t)$. This would only require storing $K$ primal solutions, one for each time segment. Transcription LSS requires $u$ at all $m$ time steps discretizing the time horizon $T_{1}-T_{0}$. Therefore, MSS will require less memory than transcription LSS since $K<m$.

Finally, the time-parallel MATVEC presented in section 5.4 allows for using distributed memory as the data needed for each time segment can be stored on different CPU’s.

Overall, the MSS formulation offers many more opportunities to minimize storage requirements than transcription LSS. However, memory efficiency alone does not make an algorithm computationally efficient. The time required to complete an algorithm is even more important in many cases.

6 Examples of MSS

6.1 Dowell’s plate

To study characteristics of the MSS algorithm including speed and accuracy, Dowell’s plate model is studied. This model is a relatively simple aeroelastic test problem first explored by Dowell [52]. The set-up, shown in figure 2, is a plate under a compressive in-plane loading $N_{x}^{E}$ in supersonic crossflow. For certain cases, Dowell observed chaotic flutter of this plate [52].

Using linear piston theory, a non-linear thin plate model, and assuming relatively high Mach numbers, the following PDE is obtained [53, 54, 55]:

[TABLE]

where the variables are defined in table 1.

A set of ordinary differential equations (ODE’s) can be obtained from equation (43) using Galerkin’s method with the modal expansion

[TABLE]

which is consistent with simply supported boundary conditions, $w=\partial^{2}w/\partial x^{2}=0$ at $x=0$ and $x=a$ [52]. This results in the following non-dimensional equations

[TABLE]

where $A_{n}\equiv a_{n}/h$, $\lambda\equiv\rho_{\infty}U_{\infty}^{3}a^{3}/MD$, $\mu\equiv\rho_{\infty}a/m$, $R_{x}\equiv N_{x}^{E}a^{2}/D$, $P\equiv\Delta pa^{4}/Dh$, $\tau\equiv t(D/ma^{4})^{1/2}$, and a prime denotes $\partial(\ )/\partial\tau$. Also define $W\equiv w/h$ for later use. All results presented in this paper use four modes ($n=4$), which is said to be sufficient by Dowell [52].

The two parameters studied are the flow velocity parameter $\lambda$ and the load parameter $R_{x}$. As $\lambda$ is increased in the absence of external compression ($R_{x}=N_{x}=0$), the plate will flutter. On the other hand, as $R_{x}$ is increased the plate will buckle. Chaotic behavior is observed with the right combination of $\lambda$ and $R_{x}$, and occurs due to coupling in the buckling and flutter instabilities [52]. MSS is tested in this chaotic region of parameter space ($\lambda,R_{x}$).

The objective function used for this study is the variance of the dimensionless transverse deflection $W(0.75a,t)$. Also, all calculations presented for Dowell’s plate use the adjoint version of MSS.

Time integration of the governing and tangent equations was done with a 3rd order explicit Runge-Kutta scheme based on the scheme in chapter 6 of [56],

[TABLE]

where subscript $j$ corresponds to the $j$th discrete time step, $\Delta t$ is the time step size and

[TABLE]

The corresponding dual consistent adjoint time stepping scheme was used for adjoint solver for the reasons discussed in section 5.2. The time step size $\Delta t$ was selected to ensure stability and accuracy of the scheme, resulting in fairly small time steps for Dowell’s plate, which is modeled by stiff ODE’s.

Figures 3 and 5 show that MSS can compute accurate sensitivities for wide ranges of $R_{x}$. Figure 5 shows that MSS is very accurate for $R_{x}<-76.0$, fairly accurate for $-76.0<R_{x}<-68.0$ and $-60.0<R_{x}<-50.0$, and inaccurate for $-68.0<R_{x}<-60.0$. The regions where MSS is most inaccurate correspond to values of $R_{x}$ for which the attractor geometry is rapidly changing with respect to $R_{x}$, as shown in figure 5. The rapid drops in attractor dimension around $R_{x}=-65.0$ and $R_{x}=-61.0$ are in the same region where the accuracy of MSS is degraded. This is consistent with the fact that sensitivity analysis is often ill-posed when attractor topology changes rapidly333Attractor topology changes such as bifurcations are inherently non-linear phenomena, so using using the tangent or adjoint equations to compute perturbations near them is a poor approximation, since the tangent and adjoint are linear. Also, there is no guarantee time-averaged objective functions will vary smoothly near topology changes [17].

As MSS is accurate for a wide range of parameters for Dowell’s plate, this test case is used to study the properties of MSS, including the conditioning of the KKT Schur complement matrix in equation (36).

6.2 Kuramoto-Sivashinsky Equation

The second example is the Kuramoto-Sivashinsky (K-S) equation, a case previously studied with LSS [12]. The K-S equation is a 4th order, chaotic PDE, that can be used to model a number of physical phenomena [59]. Kuramoto derived the equation for angular-phase turbulence for a system of reaction-diffusion equations modeling the Belouzov-Zabotinskii reaction in three spatial dimensions [60, 61]. Sivashinsky also derived the equation to model the evolution of instabilities in a distributed plane flame front [62, 63]. In addition the K-S equation has also been shown to be a model of Poiseuille flow of a film layer on an inclined plane [64]. Numerical studies typically use the 1D version of the K-S equation:

[TABLE]

The $c$ term is added to make the system ergodic [12]. The equation was solved on the domain $0\leq x\leq 128$, with the boundary conditions:

[TABLE]

These boundary conditions make the K-S equation ergodic, unlike the periodic boundary conditions used in many studies. A typical solution of the K-S equation from a randomized initial condition is shown in figure 6. The initial condition is randomized as follows: since equation (46) is discretized with a 2nd order central difference scheme, $u(x,0)$ at each node is set to a random number drawn from a uniform distribution between -0.5 and 0.5. Time stepping was conducted with the 3rd order Runge-Kutta scheme from equation (45). A time step size of 0.1 units was used for the results presented in this paper.

The objective function $\bar{J}$ considered for the K-S equation in this paper is $u^{2}$ averaged over both space and time:

[TABLE]

Figure 7 shows sensitivities of computed by MSS. As for Dowell’s plate, MSS computes accurate sensitivities for the K-S equation, just like LSS [12].

7 Rate of Convergence of MSS

One of the main takeaways from the study of LSS for chaotic vortex shedding is that the transcription LSS implementation took a long time to compute sensitivities [14]. This was mainly due to the large number of GMRES iterations (or search directions) required to converge the KKT Schur complement enough to obtain a sufficiently accurate sensitivity. To make LSS practical for large systems like those encountered in CFD simulations, convergence needs to be faster.

The convergence rates of a wide range of iterative solvers depends on the condition number $\kappa$ of the matrix. For a symmetric positive definite matrix like the MSS Schur complement, the condition number $\kappa$ is defined as the ratio between the largest and smallest eigenvalues, $\mu_{max}$ and $\mu_{min}$. Typically, $\kappa$ is quite large for the KKT Schur complement matrices derived for both transcription LSS and MSS. Figure 10 shows a typical spectrum for the MSS KKT Schur complement with $\epsilon=0$. This matrix has a condition number $\kappa\sim 10^{14}$. Large values of $\kappa$ are problematic, because iterative solvers typically converge faster for systems with smaller $\kappa$. This means that it is desirable to have $\mu_{max}$ and $\mu_{min}$ as close in value as possible. These eigenvalues depend on the chaotic dynamics of the system being analyzed, the number of time segments $K$, the time segment lengths, $\Delta T_{i}$, and the filtering parameter $\epsilon$. The details of the relationship of $\mu_{max}$ and $\mu_{min}$ to these properties and parameters is presented using some analysis of the MSS KKT Schur complement and some results for MSS applied to Dowell’s plate.

7.1 The Largest Eigenvalue

First, the largest eigenvalue, $\mu_{max}$, is approximately (see F)

[TABLE]

assuming $\max_{i}e^{\tilde{\Lambda}^{max}_{i}\Delta T_{i}}>>1$, where $\Delta T_{i}$ is the length of time segment $i$ and the quantity $\tilde{\Lambda}^{max}_{i}$ is the largest finite time Lyapunov exponent (FTLE) in time segment $i$. It defined as a finite time approximation of the largest Lyapunov exponent and

[TABLE]

The presence of $\tilde{\Lambda}^{max}_{i}$ in (48) shows that $\mu_{max}$ is linked to the dynamics of the governing equations being studied. It also explains why the magnitude of $\mu_{max}$ is typically large. For many chaotic dynamical systems, including chaotic fluid flows, it has been observed that FTLE’s can vary widely on a strange attractor [65]. This means that in some time segments $\tilde{\Lambda}^{max}_{i}$ will be much larger than $\Lambda^{max}$, in others, it may even be negative. The former case can lead to large values of $\mu_{max}$ even for relatively small time segment lengths $\Delta T_{i}$.

Figure 10 shows that the estimated value of $\mu_{max}$ according to equation (48) matches the true maximum eigenvalue very well. Also, it confirms that $\mu_{max}$ increases as the time segment size is increased. This growth in $\mu_{max}$ with time segment length $\Delta T_{i}$ is non-uniform because of the variation of the FTLE’s on the regions of strange attractor encompassed by the MSS time horizon. Since the same solution $u(t)$ and time horizon is considered for all cases in figure 10, different choices of time segments will compute different FTLE’s. The time segment lengths are uniform, so $\Delta T_{i}=(t_{K}-t_{0})/K$. Therefore, decreasing $K$ could decrease $\mu_{max}$ in some cases if the maximum FTLE $\tilde{\Lambda}^{max}_{i}$ decreases. However, these decreases are only temporary. As $K$ is decreased (and $\Delta T_{i}$ is increased) $\tilde{\Lambda}^{max}_{i}$ will approach $\Lambda^{max}$ and $\mu_{max}$ will grow exponentially with $\tilde{\Lambda}^{max}_{i}\Delta T_{i}$.

Only uniform time segment lengths are considered in this thesis, but the variation of $\tilde{\Lambda}^{max}_{i}$ and its impact on $\mu_{max}$ may justify using non-uniform time segment lengths for MSS. Non-uniform time segment lengths could be used to control $\mu_{max}$ and perhaps even the condition number $\kappa$, depending on how non-uniform time segments impact the smallest eigenvalue $\mu_{min}$.

7.2 The Smallest Eigenvalue

When LSS is conducted for a long enough time horizon, the KKT system for both MSS and transcription LSS can become nearly rank deficient, or $\mu_{min}\approx 0$.444This is the case for MSS when the filtering parameter $\epsilon$ is 0. This issue was first observed for early applications of shadowing for noise control [66], but can arise for LSS as well, as shown by the very low values of $\mu_{min}$ in figure 10. Also the slow convergence of the KKT Schur complement system for the chaotic vortex shedding case studied by Blonigan et al. [14] suggests the Schur complement is nearly rank deficient.

The small magnitude of $\mu_{min}$ when $\epsilon=0$ is a symptom of applying LSS or any shadowing approach to a system with a quasi-hyperbolic attractor as opposed to a uniformly hyperbolic attractor. In fact, it has been shown that $\mu_{min}$ has a lower bound for uniformly hyperbolic attractors [11]. As discussed in section 2, quasi-hyperbolic attractors include some points where some of the Lyapunov covariant vectors become parallel, called homoclinic tangencies [66]. When a trajectory $u(t)$ passes very close to a homoclinic tangency, it is difficult to find corresponding shadow trajectories and therefore, shadowing directions $v(t)$ [66].

To understand the difficulty of finding shadowing directions near homoclinic tangencies, consider the $i$th checkpoint for tangent MSS, where $v(t_{i})=\textbf{v}_{i}$. If this checkpoint is on a hyperbolic region of the attractor, then $\textbf{v}_{i}$ and $v(t),t\in[t_{i},t_{i+1}]$ can be expressed as a linear combination of Lyapunov covariant vectors, which act as a basis for all $n$ dimensions of phase space at the point $u(t)$. Define the minimum angle between any two Lyapunov covariant vectors at checkpoint $i$ as $\theta_{i}$. In the limit $\theta_{i}\to 0$, checkpoint $i$ coincides with a homoclinic tangency and two of the covariant vectors are parallel, so now $n-1$ of the $n$ covariant vectors are linearly independent. Therefore $\textbf{v}_{i}$ (and therefore $v(t),t\in[t_{i},t_{i+1}]$) can only be expressed a linear combination of $n-1$ covariant vectors. This means there are only $n-1$ degrees of freedom to solve the checkpoint continuity constraint equation (29). Since this is one of the equations that makes up the KKT system, the KKT system is also rank deficient. Therefore small $\theta_{i}$, or equivalently the close proximity of a checkpoint to a homoclinic tangency, implies a low condition number $\kappa$ for the entire system.

This is problematic because longer time horizons $t_{K}-t_{0}$ will increase the likelihood of a checkpoint being close to a homoclinic tangency. This occurs because as a time horizon is increased in length, a solution $u(t)$ will pass near more homoclinic tangencies or pass closer by previously visited homoclinic tangencies. This trend was observed for a number of chaotic systems, including Dowell’s plate, which is shown in figure 12. The minimum eigenvalue of the KKT system for Dowell’s plate decreases by orders of magnitude as $K$ is increased. On the other hand, the minimum eigenvalue for the uniformly hyperbolic solenoid map shown in figure 12, stays approximately constant as $K$ is increased. This behavior is expected since hyperbolic attractors do not have homoclinic tangencies.

Finally, it should be noted that the decrease in the minimum eigenvalue due to homoclinic tangencies also explains the poor conditioning of the transcription LSS KKT system. This is because of the equivalence of MSS for $\epsilon=0$ and transcription LSS as $K\to\infty$ shown in section 4 and appendix B.

7.3 Controlling the Condition Number with the Filtering Parameter

The MSS filtering parameter, $\epsilon$, can be used to reduce the condition number of the MSS KKT system considerably. It does this by applying a constant shift $\epsilon$ to the entire eigenvalue spectrum. Define $\mu_{max,0}$ and $\mu_{min,0}$ as the maximum and minimum eigenvalues of the MSS KKT Schur complement with $\epsilon=0$, or $\uuline{B}\uuline{B}^{T}$. Then the condition number of the MSS KKT Schur complement with non-zero $\epsilon$ is:

[TABLE]

Given that $\mu_{min}$ is typically much smaller than 1 for the MSS KKT Schur complement (see figure 10), a very small value of $\epsilon$ can improve the condition number $\kappa$ by orders of magnitude. This means that non-zero $\epsilon$ can reduce the number of iterations and therefore the amount of time required by an iterative solver to solve equation (36).

This can be demonstrated on Dowell’s plate. As the MSS system is shown to be nearly indefinite in figure 10, it is solved with the Krylov subspace solver MINRES [45]. Figure 13 shows that increasing $\epsilon$ results in faster convergence of the residual ${\bf r}=\uuline{{\bf A}}\underline{\hat{\textbf{w}}}+\uuline{{\bf B}}\underline{\bf g}$. At the same time,figures (15) and (15) show that the choice of $\epsilon$ affects the sensitivities computed by MSS. These figures show that there is an optimal value of $\epsilon$ between the minimum and maximum eigenvalues $\mu_{min,0}$ and $\mu_{max,0}$.

Figure (15) shows that the optimal value of $\epsilon$ is robust and a wide range of $\epsilon$ produces results with comparable accuracy for Dowell’s plate. The errors when $\epsilon\approx\mu_{min,0}$ arise because of the poor conditioning of the minimization problem for Dowell’s plate. These errors are not observed for the K-S equation, as shown in figure 16. The errors at large values of $\epsilon$ observed in figures 15 and 15 arise because adjoint MSS becomes an ensemble adjoint approach as $\epsilon\to\infty$. For large values of $\epsilon$, consider equation (36)

[TABLE]

When $\epsilon$ becomes very large relative to the entries of $\uuline{{\bf B}}\uuline{{\bf B}}^{T}$ and $\uuline{{\bf B}}\underline{\bf g}$, equation (36) becomes

[TABLE]

This corresponds to solving the adjoint in each time segment with a terminal condition of 0, as in the ensemble adjoint method [19]. The only difference is the presence of the projection operator $P_{t}$ in MSS, but this appears to have a small effect in practice, as seen in figure 16.

The ensemble adjoint method requires a time segment around the size of the largest time scale of the system being analyzed [19]. This time scale is proportional to the largest inverse Lyapunov exponent of the system.

[TABLE]

where $\Lambda^{|min|}$ is the Lyapunov exponent with the smallest non-zero magnitude. In the case of the modified K-S equation with $c=0.0$, this is roughly $1/0.004=250$ time units [12]. Since the time segments used for the results in figure 16 are only 102.4 units, MSS suffers from the same bias as the ensemble adjoint. However, lower values of $\epsilon$ give sensitivities similar to those computed in a previous study [12].

8 Conclusion

Unlike conventional sensitivity analysis methods, LSS is able to compute useful sensitivities of long-time-averaged quantities in chaotic systems. Although other sensitivity analysis methods for chaotic systems have been developed in other communities, these are not viable approaches for high-fidelity, scale-resolving turbulence simulations like LES or DNS. The issues with these other methods include the slow convergence of the ensemble adjoint approach with the number of time horizons used, and the difficulties of extending Fokker-Planck approaches to chaotic systems with many DoFs. LSS avoids many of these issues, but the original implementation of LSS, transcription LSS, suffers from large computational costs.

There are many ways the cost of LSS can be reduced by using multiple shooting rather than a transcription approach. The MSS implementation is definitely more memory efficient than transcription LSS, and can be made to converge faster with the right choice of time checkpoints and filtering parameter. Because of this, larger chaotic dynamical systems like scale-resolving fluid flow simulations should be analyzed with an MSS implementation.

There is still work to be done on the MSS approach itself. This includes exploring other solver approaches such as MGRIT and a detailed investigation of preconditioners for MSS.

Overall, this work is imperative to realizing simulation based design for aerospace vehicles and other engineering devices which require high fidelity flow simulations for accurate analysis. Simulation based design with LSS would be very valuable, leading to gains in performance and providing a deeper understanding of the physics encountered by a wide range of vehicles and devices, from turbomachinery components to heavy-lift launch vehicles.

Appendix A Computing the Maximum Lyapunov Exponent

$\Lambda^{max}$ can be computed by the following approach, a simplification of the algorithm of [41]:

Compute $\Lambda^{max}$

Inputs: Initial condition for the governing equations $u_{0}$, Spin-up time $t_{0}$, Specified time horizon and evenly spaced checkpoints $t_{0},t_{1},...,t_{K}$ with $t_{i}-t_{i-1}=\Delta t$;

Ouputs: Maximum Lyapunov exponent $\Lambda^{max}$

Time integrate the governing equations (1) to compute $u(t)$ for the specified time horizon. 2. 2.

Set $i=1$ and set $\textbf{v}_{0}$ equal to some arbitrary unitary vector. 3. 3.

Compute $v(t_{i}^{-})$ by integrating $\frac{dv}{dt}=\frac{\partial{f}}{\partial{u}}v\;t\in(t_{i-1},t_{i}]$ with the initial condition $v(t_{i-1})=\textbf{v}_{i-1}$. 4. 4.

Compute and save $R_{i}=\frac{\|v(t_{i}^{-})\|_{2}}{\|v(t_{i-1})\|_{2}}$. 5. 5.

Set $\textbf{v}_{i}=v(t_{i}^{-})$. 6. 6.

Set $i=i+1$ and repeat steps 3 to 5 until $i=K$ 7. 7.

Solve

[TABLE]

This approximation becomes more accurate as $K\to\infty$ for any choice of $\Delta t$. [41]

Appendix B Relationship between MSS and Transcription LSS

For $\alpha=0$, the LSS minimization problem (6) becomes

[TABLE]

with the KKT equations

[TABLE]

First, for uniformly spaced checkpoints $t_{i}$ and $\epsilon=0$, the MSS minimization statement multiplied by $\Delta t=t_{i+1}-t_{i}$ is

[TABLE]

For $t_{0}=T_{0}$ and $t_{K}=T_{1}$ this is equivalent to equation (50). The limit in equation (54) is true because a Riemann sum approaches the Riemann integral as smaller partitions $\Delta t$ are considered. This also ensures that the second term on the left of (54) disappears as $\lim_{K\to\infty}\delta=0$.

The transcription LSS KKT equations are also consistent with the MSS constraint (29). The KKT equations (52) and (53) is equivalent to equations (12) and (13) when $\epsilon=0$ and continuity is enforced at checkpoints $t_{i}$.

Equations (51) and (52) imply equation (13) holds when $\alpha=0$. This is because $\langle f(u(t);s),v(t)\rangle=0$ must also hold to ensure that $dw/dt$ and therefore $w(t)$ has no component parallel to $f(u(t);s)$ and equation (52) is satisfied.

Together, all of this shows that the MSS minimization problem is equivalent to the transcription LSS minimization problem as $K\to\infty$ if the checkpoints are uniformly spaced, $\epsilon=0$ and $\alpha=0$.

Appendix C MSS Gradient

The expression to compute the gradient from $v(t)$ and $\eta(t)$ is

[TABLE]

This expression can be written as a function of $v^{\prime}(t)$ from equation (17) by using the following expression

[TABLE]

where $J_{i+1}=J(u(t_{i+1});s)$. Note that

[TABLE]

where $f_{i+1}=f(u(t_{i+1});s)$, Therefore

[TABLE]

Equation (57) is substituted into equation (55) to obtain

[TABLE]

Finally,

[TABLE]

where $T\equiv t_{K}-t_{0}$.

Appendix D Adjoint MSS Derivation

For a dynamical system with a finite number of states $n$, equation (27) can be rewritten using the notation of equation (34)

[TABLE]

where $\underline{\bf g}^{T}=({\bf g}_{1}^{T},{\bf g}_{2}^{T},...{\bf g}_{K}^{T},0)$ is a $1\times(K+1)n$ vector and

[TABLE]

Where $f_{i}\equiv f(u(t_{i});s)$ and $J_{i}=J(u(t_{i});s)$. The second term in equation (58) is defined as

[TABLE]

The tangent equation (34) and the sensitivity equation (58) can be combined as follows

[TABLE]

where $\underline{\hat{\textbf{v}}}$ and $\underline{\hat{\textbf{w}}}$ are defined as the adjoint variables. Next, rearrange equation (62) as follows

[TABLE]

One can choose $\underline{\hat{\textbf{v}}}$ and $\underline{\hat{\textbf{w}}}$ to satisfy the following adjoint equation

[TABLE]

with Schur complement:

[TABLE]

If $\underline{\hat{\textbf{v}}}$ and $\underline{\hat{\textbf{w}}}$ satisfy equation (64), then the sensitivities can be computed as follows:

[TABLE]

Using equations (21) and (61), the definitions of $\underline{\bf b}$ and ${\bf h}$, respectively, equation (66) can be rewritten as:

[TABLE]

From (64):

[TABLE]

For equations (68) and (59) to be true, the adjoint solution $\hat{w}(t)$ in time segment $i-1$ should be

[TABLE]

With this choice of $\hat{w}(t)$, equation (67) simplifies to:

[TABLE]

Sensitivities of $\bar{J}$ with respect to many parameters $s$ can be computed with a single solution of $\hat{w}(t)$.

Appendix E Tangent MSS Algorithm

Since the adjoint MSS Schur complement (36) is identical to the tangent MSS Schur complement (35) the adjoint and tangent algorithms are very similar:

Tangent MSS Solver

Inputs: Initial condition for the governing equations $u_{0}$, Spin-up time $t_{0}$, Specified time horizon and checkpoints $t_{0},t_{1},...,t_{K}$, Initial guess for $Kn\times 1$ vector $\underline{\bf w}$, which contains the Lagrange multipliers at checkpoints $1$ to $K$ (default value 0);

Ouputs: Sensitivity $d\bar{J}/ds$

Calls: MATVEC algorithm that computes $\underline{\bf R}=(\uuline{{\bf B}}\uuline{{\bf B}}^{T}+\epsilon\uuline{I})\underline{\bf w}-\beta\underline{\bf b}$ where $\underline{\bf R}$ is a $Kn\times 1$ residual vector.

Time integrate the governing equations (1) to compute $u(t)$ for the specified time horizon. 2. 2.

To form the right hand side of the linear system, $\underline{\bf b}$, use the MATVEC algorithm with $\beta=-1$ and $\underline{\bf w}=0$. 3. 3.

Use some iterative algorithm to solve equation (35). To compute the left hand side $(\uuline{{\bf B}}\uuline{{\bf B}}^{T}+\epsilon\uuline{I})\underline{\bf w}$, use MATVEC with $\beta=0$. 4. 4.

Compute the sensitivity $d\bar{J}/ds$ using equation (27).

Next, a serial MATVEC algorithm is presented. Note that $t_{i}^{-}$ and $t_{i}^{+}$ refer to the time at checkpoint $i$ in time segments $i-1$ and $i$, respectively. For example, $v(t_{i}^{-})$ is the tangent solution at checkpoint $i$ computed in time segment $i-1$.

Serial Tangent MATVEC Algorithm

Inputs: $\underline{\bf w}$, a $Kn\times 1$ vector of the Lagrange multipliers at checkpoints 1 to $K$; $\beta$, a scalar;

Ouputs: $\underline{\bf R}$, a $Kn\times 1$ residual vector

MATVEC computes $\underline{\bf R}=(\uuline{{\bf B}}\uuline{{\bf B}}^{T}+\epsilon\uuline{I})\underline{\bf w}-\beta\underline{\bf b}$

For all time segments, compute $w(t_{i-1}^{+})$ by integrating $\frac{dw}{dt}=-\left(\frac{\partial{f}}{\partial{u}}\right)^{*}w,\;t\in(t_{i-1},t_{i})$ backwards in time with the terminal condition $w(t_{i})=P_{t_{i}}\textbf{w}_{i}$. 2. 2.

Save $\textbf{v}_{i-1}\equiv w(t_{i-1}^{+})-\textbf{w}_{i-1}$. If $i=1$, save $\textbf{v}_{0}\equiv\hat{w}(t_{0}^{+})$. 3. 3.

For all time segments, compute $v^{\prime}(t_{i}^{-})$ by integrating $\frac{dv^{\prime}}{dt}=\frac{\partial{f}}{\partial{u}}v^{\prime}+\beta\frac{\partial{f}}{\partial{s}},\;t\in(t_{i-1},t_{i}]$ with the initial condition $v^{\prime}(t_{i-1})=\textbf{v}_{i-1}$. 4. 4.

For all time segments, compute ${\bf R}_{i}=P_{t_{i}}\,v^{\prime}(t_{i}^{-})-\textbf{v}_{i}+\epsilon\textbf{w}_{i}$. If $i=K$, $\textbf{v}_{K}=-\textbf{w}_{K}$.

Appendix F Estimating the largest eigenvalue of the MSS KKT Schur complement

The largest eigenvalue of the MSS KKT Schur complement, $\mu_{max}$, is related to the positive Lyapunov exponent discussed in section 2.1. This can be shown with the singular value decomposition (SVD) of the tangent transition matrix $\Phi_{i}$ defined in equation (20),

[TABLE]

where ${\bf U}_{i}$ and ${\bf V}_{i}$ are orthonormal matrices and $\Sigma_{i}$ is a diagonal matrix with singular values $\sigma_{i}^{j}$ along the main diagonal. The singular values are indexed in descending order of magnitude, that is $\sigma_{i}^{1}>\sigma_{i}^{2}>\ldots>\sigma_{i}^{n}$.

Interestingly the singular value $\sigma_{i}^{1}$ can be interpreted as a finite time approximation of the largest Lyapunov exponent. The maximum Lyapunov exponent can be defined as [67]

[TABLE]

where $\sigma^{j}(\phi^{0,t})$ is the $j$th singular value of a the tangent propagator $\phi^{0,t}$, which is a $n\times n$ matrix for a system with $n$ states. Recall that $\Phi_{i}=P_{t_{i}}\phi^{t_{i-1},t_{i}}$. The projection operator $P_{t_{i}}$ only removes any component parallel to $f(u)$, which happens to be the Lyapunov covariant vector for the zero Lyapunov exponent. Therefore, $P_{t_{i}}$ has no impact on the covariant vector corresponding to $\Lambda_{max}$ and

[TABLE]

where $\Delta T_{i}=t_{i}-t_{i-1}$ is the time segment length for segment $i$. Since $\sigma_{i}^{1}$ is defined as the maximum singular value of $\Phi_{i}$ in equation (82),

[TABLE]

The quantity $\tilde{\Lambda}^{max}_{i}$ can be interpreted as a finite time approximation of $\Lambda^{max}$ in time segment $i$. From equation (84), the singular value $\sigma^{1}_{i}$ can be written in terms of the finite time Lyapunov exponent

[TABLE]

Recall from section 7.1 that $\tilde{\Lambda}^{max}_{i}\Delta t_{i}$ can be much larger than 1 due to variations in $\tilde{\Lambda}^{max}_{i}$ over each time segment. Because of this, it can be assumed that

[TABLE]

if $K$ or $\delta t_{i}$ is sufficiently large.

Next, the connection between $\mu_{max}$ and the $\sigma^{1}_{i}$ is considered. Recall from equations (35) and (26) that the MSS Schur complement matrix $\uuline{{\bf A}}$ is a block-tridiagonal matrix with $\Phi_{i}\Phi_{i}^{T}$ the main block diagonal and $\Phi_{i}$ and $\Phi_{i}^{T}$. Therefore, the eigenvalues of $\uuline{{\bf A}}$ must be related to the singular values $\sigma_{i}^{j}$ of $\Phi_{i}$. It can be shown that the largest eigenvalue $\mu_{max}$ of the $\uuline{{\bf A}}$ is

[TABLE]

To show this, consider the case where the maximum singular value is $\sigma_{l}$, corresponding the tangent propagator $\Phi_{l}$. Then, for $\underline{\bf x}=(0,\vdots,0,[{\bf U}_{l}^{1}]^{T},0,\vdots,0)^{T}$,

[TABLE]

Rows $l-1$ and $l+1$ of $\uuline{A}\underline{\bf x}$ are negligible relative to row $i$ if $\sigma_{l}^{1}\gg 1$ as in equation (86) and $\epsilon\geq 0$.

In row $l-1$, since by definition $\|{\bf U}_{l}^{1}\|_{2}=\|{\bf V}_{l}^{1}\|_{2}=1$,

[TABLE]

The magnitude of row $l+1$ of $\uuline{A}\underline{\bf x}$ is a maximum when ${\bf U}_{1}^{l}={\bf V}_{l+1}^{1}$, as ${\bf V}_{l+1}^{1}$ will be stretched by the largest singular value for time segment $l+1$, $\sigma_{l+1}^{1}$, so

[TABLE]

since $\sigma_{l}^{1}$ is the largest singular value for all time segments, $\sigma_{l}^{1}>\sigma_{l+1}^{1}$ and since by definition $\|{\bf U}_{l}^{1}\|_{2}=\|{\bf U}_{l+1}^{1}\|_{2}=1$

[TABLE]

The inequalities in equations (123) and (124) show that

[TABLE]

Therefore, the largest eigenvalue of $\uuline{A}$ is

[TABLE]

with the corresponding eigenvector $\underline{\bf x}$. Since $\sigma_{l}^{1}$ is by definition the largest singular value for all time segments, equations (126) and (87) are equivalent.

Note that $\mu_{max}\gg 1$ in the typical spectrum shown in figure 10. Since this spectrum was computed for $\epsilon=0$, the assumption $\sigma_{l}^{1}\gg 1$ is consistent with a typical spectrum of $\uuline{{\bf A}}$.

Therefore,

[TABLE]

Acknowledgments

Research for this paper was conducted with Government support under FA9550-11-C-0028 and awarded by the Department of Defense, Air Force Office of Scientific Research, National Defense Science and Engineering Graduate (NDSEG) Fellowship, 32 CFR 168a.