Optimal bounds and extremal trajectories for time averages in nonlinear dynamical systems

TL;DR

This paper establishes a duality framework for bounding long-time averages in nonlinear dynamical systems, showing auxiliary functions can provide arbitrarily sharp bounds and identify regions containing extremal trajectories.

Contribution

It proves the strong duality between extremal trajectory problems and auxiliary function optimization, enabling precise bounds and trajectory localization in nonlinear systems.

Findings

Auxiliary functions yield arbitrarily sharp upper bounds on time averages.

Nearly minimal auxiliary functions identify phase space regions containing extremal trajectories.

Semidefinite programming can construct auxiliary functions for polynomial systems like the Lorenz system.

Abstract

For any quantity of interest in a system governed by ordinary differential equations, it is natural to seek the largest (or smallest) long-time average among solution trajectories, as well as the extremal trajectories themselves. Upper bounds on time averages can be proved a priori using auxiliary functions, the optimal choice of which is a convex optimization problem. We prove that the problems of finding maximal trajectories and minimal auxiliary functions are strongly dual. Thus, auxiliary functions provide arbitrarily sharp upper bounds on time averages. Moreover, any nearly minimal auxiliary function provides phase space volumes in which all nearly maximal trajectories are guaranteed to lie. For polynomial equations, auxiliary functions can be constructed by semidefinite programming, which we illustrate using the Lorenz system.

Figures

Tables

Table 1. Table 1: Upper bounds z 4 ¯ ≤ U ¯ superscript 𝑧 4 𝑈 \overline{z^{4}}\leq U in the Lorenz system computed using polynomial V ​ ( x , y , z ) 𝑉 𝑥 𝑦 𝑧 V(x,y,z) of various degrees. Underlined digits agree with the value z 4 ¯ ≈ 592827.338 ¯ superscript 𝑧 4 592827.338 \overline{z^{4}}\approx 592827.338 attained on the shortest periodic orbit.

Degree of $V$	Upper bound $U$
4	635908.
6	595152.
8	592935.
10	592827.568
12	592827.344

Equations

$$\frac{d}{dt}\mathbf{x}=\mathbf{f}(\mathbf{x}),$$

$$\overline{\Phi }({\mathbf{x}}_0)=\underset{T\to \infty }{\mathrm{limsup}}\frac{1}{T}{\int }_0^T\Phi (\mathbf{x}(t))dt.$$

$${\overline{\Phi }}^{\ast }=\underset{{\mathbf{x}}_0\in B}{\max }\overline{\Phi }({\mathbf{x}}_0).$$

$$\overline{\frac{d}{dt}V}=\overline{\mathbf{f}\cdot \nabla V}=0.$$

$$\overline{\Phi }=\overline{\Phi +\mathbf{f}\cdot \nabla V}.$$

$$\overline{\Phi }({\mathbf{x}}_0)\le \underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}$$

$$\underset{{\mathbf{x}}_0\in B}{\max }\overline{\Phi }\le \underset{V\in C^1(B)}{\inf }\underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}.$$

$$\underset{{\mathbf{x}}_0\in B}{\max }\overline{\Phi }=\underset{V\in C^1(B)}{\inf }\underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}.$$

$$\frac{d}{dt}V=\mathbf{f}\cdot \nabla V\le U-\Phi \to -\infty$$

$$\overline{\Phi }({\mathbf{x}}_0^{\ast })=\underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V^{\ast }\right\}.$$

$$0\le \underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}-\overline{\Phi }({\mathbf{x}}_0)\le ϵ.$$

$$0\le \overline{\underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}-(\Phi +\mathbf{f}\cdot \nabla V)}\le ϵ$$

$$\mathcal{S}_{M}=\Big{\{}\mathbf{x}\in B:\max_{\mathbf{x}\in B}\left\{\Phi+\mathbf{f}\cdot\nabla V\right\}-(\Phi+\mathbf{f}\cdot\nabla V)(\mathbf{x})\leq M\Big{\}}.$$

$$\underset{T\to \infty }{\mathrm{liminf}}{\mathcal{F}}_M(T)\ge 1-ϵ/M.$$

$$ϵ\ge \overline{M\text{mathbbm}{1}_{\mathbf{x}\notin {\mathcal{S}}_M}}=M\left(1-\underset{T\to \infty }{\mathrm{liminf}}{\mathcal{F}}_M(T)\right).$$

$$\Phi (\mathbf{x})+\mathbf{f}(\mathbf{x})\cdot \nabla V^{\ast }(\mathbf{x})={\overline{\Phi }}^{\ast }$$

$$\underset{{\mathbf{x}}_0\in B}{\max }\overline{\Phi }$$

$$=\underset{\mu \in Pr(B)}{\sup }\underset{V\in C^1(B)}{\inf }\int \Phi +\mathbf{f}\cdot \nabla Vd\mu$$

$$=\underset{V\in C^1(B)}{\inf }\underset{\mu \in Pr(B)}{\sup }\int \Phi +\mathbf{f}\cdot \nabla Vd\mu$$

$$=\underset{V\in C^1(B)}{\inf }\underset{\mathbf{x}\in B}{\max }\left\{\Phi +\mathbf{f}\cdot \nabla V\right\}.$$

$$\underset{{\mathbf{x}}_0\in B}{\sup }\overline{\Phi }({\mathbf{x}}_0)\le \underset{\begin{array}{c}\mu \in \mathrm{Pr}(B)\\ \mu \text{ is invar.}\end{array}}{\max }\int \Phi d\mu .$$

$$\overline{\Phi }({\mathbf{x}}_0)=\int \Phi d{\mu }^{\ast }=\underset{\begin{array}{c}\mu \in \mathrm{Pr}(B)\\ \mu \text{ is invar.}\end{array}}{\max }\int \Phi d\mu$$

$$\int \mathbf{f}\cdot \nabla \psi d\mu =0$$

$$\frac{d}{dt}\int \psi \circ {\phi }_{t}d\mu =\int \mathbf{f}\cdot \nabla (\psi \circ {\phi }_t)d\mu$$

$$\underset{V\in C^1(B)}{\inf }\int \mathbf{f}\cdot \nabla Vd\mu =\{\begin{array}{ll}0& \text{div}\mathbf{f}\mu =0\\ -\infty & \text{otherwise}.\end{array}$$

$$\underset{\mu \in Pr(B)}{\sup }\underset{V\in C^1(B)}{\inf }\int \Phi +\mathbf{f}\cdot \nabla Vd\mu =\underset{\begin{array}{c}\mu \in \mathrm{Pr}(B)\\ \mu \text{ is invar.}\end{array}}{\max }\int \Phi d\mu .$$

$$\mathcal{L}(\mu ,V)=\int \Phi +\mathbf{f}\cdot \nabla Vd\mu$$

$$\underset{\mu \in \mathrm{Pr}(B)}{\sup }\underset{V\in C^1(B)}{\inf }\mathcal{L}=\underset{V\in C^1(B)}{\inf }\underset{\mu \in \mathrm{Pr}(B)}{\sup }\mathcal{L}.$$

$$\frac{d}{dt}\mathbf{u}=\mathbf{f}\{\mathbf{u}\}.$$

Full text

Keywords: Nonlinear dynamical systems, Time averages, Ergodic theory, Semidefinite programming, Sum-of-squares polynomials, Lorenz equations

1. Introduction  For dynamical systems governed by ordinary differential equations (ODEs) whose solutions are complicated and perhaps chaotic, the primary interest is often in long-time averages of key quantities. Time averages can depend on initial conditions, so it is natural to seek the largest or smallest average among all trajectories, as well as the extremal trajectories that realize them. For various purposes including the control of chaos Yang et al. (2000), it is valuable to know extremal trajectories regardless of their stability. In other situations one is interested only in stable trajectories, but determining extrema only among these can be prohibitively difficult. The next best option is to determine extrema among all trajectories.

One common way to seek extremal time averages is to construct a large number of candidate trajectories. However, for many nonlinear systems it is challenging both to compute trajectories and to determine that the extremal ones have not been overlooked. In this Letter we study an alternative approach that is broadly applicable and often more tractable: constructing sharp a priori bounds on long-time averages. We focus on upper bounds; lower bounds are analogous.

The search for an upper bound on a long-time average can be posed as a convex optimization problem Chernyshenko et al. (2014), as described in the next section. Its solution requires no knowledge of trajectories. What is optimized is an auxiliary function defined on phase space, similar to but distinct from Lyapunov functions in stability theory. We prove here that the best bound produced by solving this convex optimization problem coincides exactly with the extremal long-time average. That is, arbitrarily sharp bounds on time averages can be produced using increasingly optimal auxiliary functions. Moreover, nearly optimal auxiliary functions yield volumes in phase space where maximal and nearly maximal trajectories must reside. Whether such auxiliary functions can be computed in practice depends on the system being studied, but when the ODE and quantity of interest are polynomial, auxiliary functions can be constructed by solving semidefinite programs (SDPs) Chernyshenko et al. (2014); Fantuzzi et al. (2016); Goluskin (2017). The resulting bounds can be arbitrarily sharp. We illustrate these methods using the Lorenz system Lorenz (1963).

Consider a well-posed autonomous ODE on $\mathbb{R}^{d}$,

[TABLE]

whose solutions are continuously differentiable in their initial conditions. To guarantee this, we assume that $\mathbf{f}(\mathbf{x})$ is continuously differentiable. Given a continuous quantity of interest $\Phi(\mathbf{x})$, we define its long-time average along a trajectory $\mathbf{x}(t)$ with initial condition $\mathbf{x}(0)=\mathbf{x}_{0}$ by

[TABLE]

Time averages could be defined using $\liminf$ instead; our results hold mutatis mutandis 111The $\limsup$ and $\liminf$ averages need not coincide on every trajectory, but their maxima over trajectories do..

Let $B\subset\mathbb{R}^{d}$ be a closed bounded region such that trajectories beginning in $B$ remain there. In a dissipative system $B$ could be an absorbing set; in a conservative system $B$ could be defined by constraints on invariants. We are interested in the maximal long-time average among all trajectories eventually remaining in $B$:

[TABLE]

As shown below, there exist $\mathbf{x}_{0}$ attaining the maximum. The fundamental questions addressed here are: what is the value of $\overline{\Phi}^{*}$, and which trajectories attain it?

2. Bounds by convex optimization  Upper bounds on long-time averages can be deduced using the fact that time derivatives of bounded functions average to zero. Given any initial condition $\mathbf{x}_{0}$ in $B$ and any $V(\mathbf{x})$ in the class $C^{1}(B)$ of continuously differentiable functions on $B$ 222Here $C^{1}(B)$ denotes functions on $B$ admitting a continuously differentiable extension to a neighborhood of $B$.,

[TABLE]

This generates an infinite family of functions with the same time average as $\Phi$ since for all such $V$

[TABLE]

Bounding the righthand side pointwise gives

[TABLE]

for all initial conditions $\mathbf{x}_{0}\in B$ and auxiliary functions $V\in C^{1}(B)$. Expression (6) is useful since no knowledge of trajectories is needed to evaluate the righthand side.

To obtain the optimal bound implied by (6), we minimize the righthand side over $V$ and maximize the lefthand side over $\mathbf{x}_{0}$:

[TABLE]

The minimization over auxiliary functions $V$ in (7) is convex, although minimizers need not exist. The main mathematical result of this Letter is that the lefthand and righthand optimizations are dual variational problems, and moreover that strong duality holds, meaning that (7) can be improved to an equality:

[TABLE]

Thus, arbitrarily sharp bounds on the maximal time average $\overline{\Phi}^{*}$ can be obtained using increasingly optimal $V$.

The auxiliary function method is not the same as the various Lyapunov-type methods used to show stability or boundedness in ODE systems. However, in instances where $\Phi(\mathbf{x})$ approaches infinity as $|\mathbf{x}|\to\infty$, auxiliary functions that imply finite upper bounds $\overline{\Phi}\leq U$ also imply the existence of trapping sets by the following argument. Suppose $V\in C^{1}(\mathbb{R}^{d})$ is an auxiliary function for which the maximum of $\Phi+\mathbf{f}\cdot\nabla V$ over $\mathbb{R}^{d}$ is no larger than $U$. Then,

[TABLE]

as $|\mathbf{x}|\to\infty$. Expression (9) is a typical Lyapunov-type condition implying that all sufficiently large sublevel sets of $\Phi$ must be trapping sets.

The remainder of this Letter is organized as follows. The next section describes how nearly optimal $V$ can also be used to locate maximal and nearly maximal trajectories in phase space. The section after illustrates these ideas using the Lorenz system, for which we have constructed nearly optimal $V$ by solving SDPs. The final section proves the strong duality (8) and establishes the existence of maximal trajectories.

3. Near optimizers  In light of the duality (8), an initial condition $\mathbf{x}_{0}^{*}$ and auxiliary function $V^{*}$ are optimal if and only if they satisfy

[TABLE]

Even if the infimum over $V$ in (8) is not attained, there exist nearly optimal pairs. That is, for all $\epsilon>0$ there exist $(\mathbf{x}_{0},V)$ for which (6) is within $\epsilon$ of an equality:

[TABLE]

In such cases, $\max_{\mathbf{x}\in B}\left\{\Phi+\mathbf{f}\cdot\nabla V\right\}$ is within $\epsilon$ of being a sharp upper bound on $\overline{\Phi}^{*}$, while the trajectory starting at $\mathbf{x}_{0}$ achieves a time average $\overline{\Phi}$ within $\epsilon$ of $\overline{\Phi}^{*}$.

Nearly optimal $V$ can be used to locate all trajectories consistent with (11). Moving the constant term inside the time average and subtracting the identity (4) gives

[TABLE]

for such trajectories. The integrand in (12) is nonnegative, and the fraction of time it exceeds $\epsilon$ can be estimated. Consider the set where the integrand is no larger than $M>\epsilon$,

[TABLE]

Let $\mathcal{F}_{M}(T)$ denote the fraction of time $t\in[0,T]$ during which $\mathbf{x}(t)\in\mathcal{S}_{M}$. For any trajectory obeying (12), this time fraction is bounded below as

[TABLE]

This follows from an application of Markov’s inequality: as the integrand in (12) is nonnegative,

[TABLE]

In practice, it may not be known if there exist trajectories satisfying (11) for a given $V$ and $\epsilon$. Still, the estimate (14) says that any such trajectories would lie in $\mathcal{S}_{M}$ for a fraction of time no smaller than $1-\epsilon/M$. The conclusion is strongest when $\epsilon\ll M$, but if $M$ is too large the volume $\mathcal{S}_{M}$ is large and featureless, failing to distinguish nearly maximal trajectories. The result is most informative when $V$ is nearly optimal so that there exist trajectories where $\epsilon\ll M$ with $M$ not too large.

If a minimal $V^{*}$ exists, its set $\mathcal{S}_{0}$ is related to maximal trajectories. Any such trajectory achieves $\epsilon=0$ in (12). If it is a periodic orbit, for instance, it must lie in $\mathcal{S}_{0}$. Thus $V^{*}$ is determined up to a constant on maximal orbits. More generally, $V^{*}$ must satisfy

[TABLE]

for all $\mathbf{x}\in\mathcal{S}_{0}$. It is tempting to conjecture that $\mathcal{S}_{0}$ coincides with maximal trajectories but, as described at the end of the next section, $\mathcal{S}_{0}$ can also contain points not on any maximal trajectory.

4. Nearly optimal bounds and orbits in the Lorenz system  When $\mathbf{f}(\mathbf{x})$ and $\Phi(\mathbf{x})$ are polynomials, $V(\mathbf{x})$ can be optimized computationally within a chosen polynomial ansatz by solving an SDP Chernyshenko et al. (2014); Fantuzzi et al. (2016); Goluskin (2017). The bound $\overline{\Phi}^{*}\leq U$ follows from (6) if $\Phi+\mathbf{f}\cdot\nabla V\leq U$ for all $\mathbf{x}\in B$. A sufficient condition for this is that $U-\Phi-\mathbf{f}\cdot\nabla V$ is a sum of squares (SOS) of polynomials. The latter is equivalent to an SDP and is often computationally tractable Parrilo (2003); Lasserre (2010).

It does not follow from the strong duality result (8) that bounds computed by SOS methods can be arbitrarily sharp. This is because requiring the polynomial $U-\Phi-\mathbf{f}\cdot\nabla V$ to be SOS is generally stronger than requiring it to be nonnegative Hilbert (1888); Parrilo (2003). Nonetheless, in the few examples where time averages have been bounded using SOS methods Fantuzzi et al. (2016); Goluskin (2017), the bounds either are sharp or appear to become sharp as the polynomial degree of $V$ increases.

The remainder of this section presents the results of SOS bounding computations for the Lorenz system at the standard chaotic parameters $(\beta,\sigma,r)=(8/3,10,28)$. We obtain nearly sharp bounds on the maximal time average of $\Phi(x,y,z)=z^{4}$, as well as approximations to maximal trajectories. Because there exist compact absorbing balls Lorenz (1979), maximization over such $B$ in (3) is equivalent to maximization over $\mathbb{R}^{d}$. As reported in Goluskin (2017), searching among the periodic orbits computed by Viswanath Viswanath (2003) suggests that the maximal average $\overline{z^{4}}^{*}$ is attained by the shortest periodic orbit—the black curves in Fig. 1. We have used SOS methods to construct nearly optimal $V(x,y,z)$ and accompanying upper bounds $U$. Similar results for various $\Phi$ in the Lorenz system appear in Goluskin (2017), along with a more detailed discussion of computational implementation. Here we report more precise computations for $\Phi=z^{4}$, obtained using the multiple precision SDP solver SDPA-GMP Yamashita et al. (2010); Nakata (2010). Conversion of SOS conditions to SDPs was automated by YALIMP Löfberg (2004, 2009), which was interfaced with the solver via mpYALMIP Fantuzzi (2017).

Table 1 reports upper bounds computed by solving SDPs that produce optimal $V$ of various polynomial degrees. As the degree of $V$ increases, the bounds approach the value of $\overline{z^{4}}$ on the shortest periodic orbit to within 7 significant figures. This suggests that $\overline{z^{4}}$ is maximized on this orbit, and it reflects the sharpness of the bounds asserted by the duality (8). We do not report the lengthy expressions for these $V$; some simpler examples appear in Goluskin (2017).

To demonstrate how the volumes $\mathcal{S}_{M}$ defined in (13) approximate maximal trajectories, we consider the polynomials $V$ of degrees 6 and 10 that produce the bounds in Table 1. For the maximum in the definition of $\mathcal{S}_{M}$ we use the corresponding $U$, which bounds it from above. In each case we find that $U$ is within $0.1$ of the true maximum over any ball $B$ containing the attractor.

Figure 1a shows the volume $\mathcal{S}_{3000}$ for the degree-6 $V$, as well as the orbit that appears to maximize $\overline{z^{4}}$. The volume captures the rough location and shape of the orbit while omitting much of the strange attractor, but this $V$ is not optimal enough to yield strong quantitative statements. It follows from (14) that any trajectory where $\overline{z^{4}}$ is within $\epsilon$ of the upper bound $U=595152$ must lie inside $\mathcal{S}_{3000}$ for a fraction of time no less than $1-\epsilon/3000$. However, there are no trajectories on which this is close to unity; the higher-degree bounds in Table 1 preclude any trajectories with $\epsilon<2324$.

The degree-10 $V$ gives a significantly refined picture of maximal and nearly maximal trajectories for $\overline{z^{4}}$. Figure 1b shows the volume $\mathcal{S}_{1000}$ defined using this $V$. It follows from (14) that any trajectory where $\overline{z^{4}}$ comes within $\epsilon$ of $U=592827.568$ must lie in $\mathcal{S}_{1000}$ for a fraction of time no less than $1-\epsilon/1000$. There exist trajectories on which this is nearly unity: on the shortest periodic orbit $\overline{z^{4}}$ is only $\epsilon\approx 0.23$ smaller than $U$. Any trajectory where $\overline{z^{4}}$ is so large must spend at least 99.97% of its time in $\mathcal{S}_{1000}$.

Finding maximal trajectories directly may be intractable in many systems. We propose that the next best option is to compute volumes like those in Fig. 1. However, we caution that finding points in a set $\mathcal{S}_{M}$ defined by (13) can itself be difficult, even for polynomials.

As the auxiliary functions producing upper bounds on $\overline{\Phi}^{*}$ approach optimality, the integrand in (12) approaches zero almost everywhere on maximal trajectories. This can be seen in Fig. 2a, where the integrand is plotted along the shortest periodic orbit in the Lorenz system for our polynomials $V$ of degrees 6, 8, and 10. Along other orbits where $\overline{z^{4}}$ is large but not maximal, $V$ is less strongly constrained. As an example, we consider the periodic orbit computed in Viswanath (2003) that winds around the two wings of the Lorenz attractor with symbol sequence $AABABB$. On this orbit $\overline{z^{4}}$ is smaller than the maximum by approximately 2798. The integral in (12) remains between 0 and 2798 as $V$ approaches optimality but need not approach 0 on this orbit. In our computations it does not, as seen in Fig. 2b.

Although the auxiliary polynomials $V$ yielding the bounds on $\overline{z^{4}}$ in Table 1 approach optimality, they are not exactly optimal. Optimal $V^{*}$ which are polynomial have been constructed to prove sharp bounds on other averages in the Lorenz system, including $\overline{z}$, $\overline{z^{2}}$, and $\overline{z^{3}}$ Malkus (1972); Knobloch (1979); Goluskin (2017). These averages are maximized on the two nonzero equilibria; in each case the set $\mathcal{S}_{0}$ corresponding to $V^{*}$ is the line through these equilibria. These $\mathcal{S}_{0}$ notably include points not on any maximal trajectory. In contrast, for $\overline{z^{4}}$ the shortest periodic orbit appears to be maximal. This conjecture could be proved by constructing a $V$ whose $\mathcal{S}_{0}$ contains the shortest orbit. If such a $V$ exists, it would necessarily be optimal. However, we expect that this orbit is non-algebraic and that no polynomial $V$ can be optimal.

5. Proof of duality  To prove the strong duality (8) we require several facts from ergodic theory, which are provable by standard methods as in Walters (1982). (See also (Phelps, 2001, Chap. 12).) Let $\varphi_{t}(\mathbf{x})$ denote the flow map $\mathbf{x}(\cdot)\mapsto\mathbf{x}(\cdot+t)$ for the ODE (1). By assumption, $\varphi_{t}$ is well-defined on $B$ for all $t\geq 0$ and is continuously differentiable there. Let $Pr(B)$ denote the space of Borel probability measures on $B$. A measure $\mu\in Pr(B)$ is invariant with respect to $\varphi_{t}$ if $\mu(\varphi_{t}^{-1}A)=\mu(A)$ for all Borel sets $A$ and all $t$. Such a measure is ergodic if to any invariant Borel set it assigns measure either zero or one. The set of invariant probability measures on $B$ is nonempty, convex, and weak-$*$ compact; its extreme points are ergodic.

Our proof of the duality (8) proceeds via a standard minimax template from convex analysis (see, e.g., Ekeland and Témam (1999)). It suffices to establish the following sequence of equalities:

[TABLE]

In (17a) we reformulate our problem as a maximization over invariant measures, whose analogue for discrete maps is the topic of the field of ergodic optimization Jenkinson (2006). The remainder of this section is devoted to proving the first three equalities (17a)–(17c), along with the fact that the maximum in (17a) is attained. The final equality (17d) is evident since, for each $V$, the supremum in (17c) is attained by a suitable Dirac measure.

We begin by proving (17a). We claim that the right hand problem appearing there is a concave relaxation of the lefthand problem, and that it attains the same maximum. To see this, note first that for each initial condition $\mathbf{x}_{0}$ in $B$ there exists an invariant probability measure $\mu$ that attains $\overline{\Phi}(\mathbf{x}_{0})=\int\Phi\,d\mu$. Thus,

[TABLE]

The righthand problem in (18) is a maximization of a continuous linear functional over a compact convex subset of $Pr(B)$, so it achieves its maximum at an extreme point $\mu^{*}$ (Lax, 2002, Chap. 13), which is an ergodic invariant measure. By Birkhoff’s ergodic theorem Walters (1982),

[TABLE]

for almost every $\mathbf{x}_{0}$ in the support of $\mu^{*}$. Therefore the inequality in (18) is in fact an equality, and any such $\mathbf{x}_{0}$ attains the maximal time average $\overline{\Phi}^{*}$. This proves (17a).

To prove the second equality (17b) we require the following equivalence of Lagrangian and Eulerian notions of invariance: a Borel probability measure $\mu$ is invariant with respect to $\varphi_{t}$ by the usual (Lagrangian) definition if and only if the vector-valued measure $\mathbf{f}\mu$ is weakly divergence-free. The latter condition, which we denote by $\text{div}\,\mathbf{f}\mu=0$, means that

[TABLE]

for all smooth and compactly supported $\psi(\mathbf{x})$. This is an Eulerian characterization of invariance.

The fact that $\text{div}\,\mathbf{f}\mu=0$ is equivalent to invariance is quickly proved using the flow semigroup identity, which states that $\varphi_{t+s}=\varphi_{t}\circ\varphi_{s}$ for all $t$ and $s$. It follows that

[TABLE]

for all smooth and compactly supported $\psi$. If $\text{div}\,\mathbf{f}\mu=0$, the righthand side of (21) vanishes, so $\mu$ is invariant. Conversely, if $\mu$ is invariant then the lefthand side of (21) vanishes for all $t$, and at $t=0$ we find the statement that $\mathbf{f}\mu$ is weakly divergence-free.

With the Eulerian characterization of invariance in hand, we turn to proving (17b). Depending on $\mu$, there are two possibilities for the minimization over $V$ in (17b):

[TABLE]

Only measures for which $\text{div}\,\mathbf{f}\mu=0$ can give values larger than $-\infty$ in (17b). As shown above, $\text{div}\,\mathbf{f}\mu=0$ if and only if $\mu$ is invariant. Therefore, since there always exists at least one invariant probability measure,

[TABLE]

Thus (17b) is proven. In other words,

[TABLE]

is a Lagrangian for the constrained maximization appearing on the righthand side of (23).

Finally, we prove the equality (17c). In terms of the Lagrangian $\mathcal{L}$, we must show that

[TABLE]

The fact that the order of $\inf$ and $\sup$ can be reversed without introducing a so-called duality gap is not trivial; it is at the heart of our proof of the strong duality (8). This reversal relies on properties of the Lagrangian $\mathcal{L}$ and the spaces $Pr(B)$ and $C^{1}(B)$.

The desired equality (25) can be proved using any of several abstract minimax theorems from convex analysis. Here we apply a fairly general infinite-dimensional version due to Sion Sion (1958). We follow the notation of its statement in the introduction of Komiya (1988), which contains an elementary proof. Let $X=\Pr(B)$ in the weak-$*$ topology. It is a compact convex subset of a linear topological space. Let $Y=C^{1}(B)$ in the $C^{1}$-norm topology, which is itself a linear topological space. Take $f=-\mathcal{L}$ and observe that $f(x,\cdot)$ is upper semicontinuous and quasi-concave on $Y$ for each $x\in X$, and that $f(\cdot,y)$ is lower semicontinuous and quasi-convex on $X$ for each $y\in Y$. Then (25) follows from a direct application of Sion’s minimax theorem Komiya (1988), so (17c) is proven.

This completes the proof of the equalities (17a)–(17d) and so too the proof of the strong duality (8).

6. Conclusions  This Letter establishes that the auxiliary function method for proving a priori bounds on long-time averages in dynamical systems yields arbitrarily sharp bounds, so long as the dynamics arise from ODEs. The proof elucidates the role that auxiliary functions play in the search for optimal bounds: they are Lagrange multipliers enforcing the constraint of invariance for probability measures on phase space. We also have demonstrated that certain sets constructed from nearly optimal auxiliary functions can be used to locate all optimal and nearly optimal trajectories. How close an auxiliary function is to optimality determines the fraction of time nearly optimal trajectories are guaranteed to spend in these sets. We expect much can be learned about the shape of optimal trajectories and their invariant measures by the auxiliary function approach.

Many of these observations extend to infinite-dimensional dynamics that are governed by nonlinear partial differential equations (PDEs) of the form

[TABLE]

Auxiliary functionals $V\{\mathbf{u}\}$ defined on a suitable function space yield a priori bounds on long-time averages just as in the finite-dimensional case. The “background method” used to bound mean quantities in fluid dynamics and other systems Doering and Constantin (1992) is an example of using quadratic $V$ Chernyshenko (2017). Whether or not nearly optimal functionals are always guaranteed to exist, and if they are ever quadratic for systems of interest, remains unclear. This emphasizes the need for a rigorous proof of duality between auxiliary functionals and extremal trajectories for general PDEs.

Acknowledgements  We thank Lora Billings, Rich Kerswell, Edward Ott, Ralf Spatzier, Divakar Viswanath, and Lai-Sang Young for helpful discussions and encouragement. This work was supported by NSF Award DMS-1515161, Van Loo Postdoctoral Fellowships (IT, DG), and a Guggenheim Foundation Fellowship (CRD).