Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization

We describe methods for proving upper and lower bounds on infinite-time averages in deterministic dynamical systems and on stationary expectations in stochastic systems. The dynamics and the quantities to be bounded are assumed to be polynomial functions of the state variables. The methods are computer-assisted, using sum-of-squares polynomials to formulate sufficient conditions that can be checked by semidefinite programming. In the deterministic case, we seek tight bounds that apply to particular local attractors. An obstacle to proving such bounds is that they do not hold globally; they are generally violated by trajectories starting outside the local basin of attraction. We describe two closely related ways past this obstacle: one that requires knowing a subset of the basin of attraction, and another that considers the zero-noise limit of the corresponding stochastic system. The bounding methods are illustrated using the van der Pol oscillator. We bound deterministic averages on the attracting limit cycle above and below to within 1%, which requires a lower bound that does not hold for the unstable fixed point at the origin. We obtain similarly tight upper and lower bounds on stochastic expectations for a range of noise amplitudes. Limitations of our methods for certain types of deterministic systems are discussed, along with prospects for improvement.


Introduction
In the study of dynamical systems with complicated and possibly chaotic dynamics, average quantities are often of more interest than any particular solution trajectory.This is partly because of the difficulty of computing a trajectory precisely and partly because average quantities are more important in many applications.For instance, one might seek time-averaged drag forces in a model of an oil pipeline or ensemble-averaged temperatures in a stochastic climate model.One way to estimate averages quantitatively is to integrate the dynamical system numerically and average over the resulting particular solution.Such direct numerical simulations are often straightforward, but the accuracy of the result is not guaranteed unless errors are rigorously controlled.Moreover, even a perfectly accurate solution does not give information about the different trajectories that result from different initial conditions.Finally, the computational cost of computing well-converged averages is often prohibitive, especially in systems that are high-dimensional or stochastic.
A complementary approach, which we pursue here, is to prove bounds on average quantities directly from the governing equations.An advantage over numerical integration is that bounds can be proven without knowing any solution trajectories.Furthermore, they can be proven for all possible trajectories at once, or for all trajectories within a given region of state space.On the other hand, it is generally difficult to prove bounds that are tight enough to give good estimates of the average quantities being bounded.
The aim of this work is to develop methods for proving bounds that are tight, meaning that the upper and lower bounds are equal or nearly so.We assume that the quantity of interest can be described by a function ϕ(x), where x(t) ∈ X is the state vector of a dynamical system, and we seek to bound averages of ϕ.In deterministic systems, we consider averages over infinite time, ϕ(x) := lim (1.1) If the above limit does not exist it can be replaced by limit superior or inferior for the upper and lower bound problems, respectively.The value of ϕ depends in general on which trajectory x(t) is being averaged over.In stochastic systems with a stationary probability distribution ρ(x), we consider stationary ensemble averages, ϕ(x) := X ϕ(x)ρ(x) dx. (1.2) One obstacle to proving tight bounds is, for reasons described shortly, the need to determine whether certain complicated expressions are sign-definite.As proposed in [8], this can be done systematically with computer assistance for finite-dimensional systems with polynomial dynamics.The main idea, described further in §2, is to construct a polynomial whose non-negativity implies the desired bound.By the methods of sum-ofsquares (SoS) programming [32,33], a sufficient condition for this non-negativity can then be posed as a semidefinite program (SDP).
In deterministic systems there is a second obstacle to proving tight bounds.It is generally easiest to construct bounds that hold for all possible initial conditions.Sometimes this is desired, but other times one is interested in a particular local attractor, and bounds holding globally are not generally tight for averages over a local attractor.In this sense, the local bound is spoiled by any other invariant structure in the state space, such as another attractor or an unstable fixed point.
We pursue two ways of obtaining tight bounds specific to a local attractor.The first is to enforce conditions that imply the bound only on an absorbing set around the attractor, thereby omitting other invariant structures.The second is to add noise to the system and prove bounds for ensemble averages in the vanishing noise limit.If the system is stochastically stable, then under certain conditions this limit will agree with the corresponding deterministic time average [42].Note that these ideas are applicable irrespectively of any special structure in the system (such as Hamiltonian), and can in principle be applied to systems of arbitrarily high but finite dimension.
Throughout this work, we illustrate the methods described using the van der Pol oscillator [22], which can be written as where a dot denotes d dt .The parameter µ > 0 sets the strength of the nonlinear damping.There is a limit cycle that attracts all trajectories except the unstable fixed point at the origin (Figure 1(a)), and the global invariant set is composed of the limit cycle and the fixed point.The system is a standard example of a nonlinear oscillator and has been studied extensively, including with stochastic forcing [24,12,43].Here we find nearly tight bounds on averages of x 2 + y 2 both with and without noise.The rest of this work is organized as follows.Section 2 reviews the framework of [8] for bounding deterministic averages and uses it to find upper bounds on x 2 + y 2 in the van der Pol system.Section 3 extends the framework to give attractor-specific bounds, which we use to find lower bounds on x 2 + y 2 over the van der Pol limit cycle.These lower bounds are larger than zero and thus do not apply to the unstable fixed point.The bounding methods for stochastic dynamical systems are described in §4, and they are specialized to the case of small and vanishing noise in § §5-6.The methods of these sections give bounds on x 2 +y 2 in the van der Pol example for a range of noise amplitudes.Section 7 discusses the limitations of our methods and gives ideas for improvement, and §8 offers concluding remarks.

Global bounds for deterministic systems
Consider an autonomous dynamical system in which all trajectories x(t) remain bounded as t → ∞.We assume nothing special about the structure of f except that the system is bounded, although eventually we restrict attention to polynomial f .We wish to prove upper and lower bounds on ϕ, the average of a function ϕ(x) over a trajectory x(t).In applications, ϕ can be chosen as a quantity of interest for the system.Unless the system has a single attractor that is globally asymptotically stable, different trajectories can give different values of ϕ.In this section we seek global bounds on ϕ, meaning they hold for every possible value of ϕ.Suppose we wish to prove a constant lower bound L on all possible ϕ, ϕ ≥ L. (2.2) Central to our method will be a suitably chosen differentiable storage function V (x).No matter what V is chosen, f • ∇V = 0 along any bounded trajectory because The desired bound (2.2) is thus equivalent to the inequality for any differentiable V (x).The above time average cannot be evaluated without knowing trajectories x(t), but a sufficient condition for the inequality to hold is for it to hold pointwise for all x.Thus, we will have proven that ϕ ≥ L if we can find any differentiable An upper bound ϕ ≤ U can be proven in a similar way by reversing the inequality sign in the pointwise sufficient condition.The following proposition summarizes both conditions.
Proposition 1.Let ẋ = f (x) with x ∈ R n be a dynamical system whose trajectories are bounded forward in time, let ϕ(x) be a scalar function and let ϕ be its time average defined as in (1.1).If there exist differentiable functions V u (x), V l (x), and constants U , L such that ) Remark 1.In the statement of Proposition 1 we have assumed that the time average ϕ exists.Should the limit in (1.1) not converge, the Proposition holds if we take the limit superior when computing the upper bound, and the limit inferior when computing the lower bound.That is, for any trajectory x(t), See [21] for a discussion of when infinite-time averages do or do not converge.
There are two difficulties in applying the above proposition to yield good bounds U and L. The first is choosing the storage functions V u and V l .The second is checking whether D u ≤ 0 and D l ≥ 0 for candidate storage functions and bound values.These difficulties can be prohibitive in general, but the task is greatly simplified when f and ϕ are polynomials of the state variables x 1 , . . ., x n .
In the rest of this work we assume that f and ϕ are polynomials.If the chosen V u and V l are also polynomials, then so are D u and D l .Checking the sufficient conditions (2.6a) and (2.6b) thus amounts to verifying the non-negativity of polynomial expressions.While this is an NP-hard problem, the computational complexity can be significantly reduced by replacing the conditions D l (x) ≥ 0 and −D u (x) ≥ 0 with the stronger conditions that D l and −D u belong to the set Σ of polynomials that are sums-of-squares (SoS).A polynomial P (x) is said to be SoS if it can be expressed as the sum of squares of some other polynomials-that is, if there exist polynomials {p i (x)} M i=1 such that To prove the bounds U and L it suffices to find polynomials V u and V l such that D l ∈ Σ and −D u ∈ Σ.The best bounds that can be proven in this framework are min (2.9) For the storage functions V u and V l , we must specify polynomial ansatze with undetermined coefficients.The decision variables in the upper bound optimization, for instance, are U and the coefficients of V u .
Computational methods.Optimization problems with SoS constraints such as (2.9) can be solved numerically using the methods of SoS programming.The main idea behind SoS programming is that every polynomial can be represented as a symmetric matrix (after defining a suitable polynomial basis set), and this matrix can be positive semidefinite if and only if the polynomial admits a SoS decomposition [32,33,4].Constraints involving SoS polynomials, including additional equality and inequality constraints, can be posed as equality and inequality conditions on symmetric matrices.An optimization problem with constraints of this type is known as a semidefinite program (SDP).A number of efficient computer solvers are available for SDPs (e.g.[35,38,39,10,3]), and the software packages YALMIP [25] or SOSTOOLS [29] can assist in formulating SoS constraints as SDPs.More details on convex optimization and the link between SoS polynomials and SDPs can be found in [5,32,26], while examples using SoS programming to study dynamical systems can be found in [32,30,31,37,14,7,16,40,18].
For our numerical implementation, we used the SoS module of YALMIP [26] to transform SoS optimization problems into SDPs.To solve the SDPs, we used the multipleprecision solver SDPA-GMP [10] for the following reasons.First, the SDPs we solved were badly conditioned even for modest polynomial degrees, and none of the standard doubleprecision solvers we tested converged reliably.This issue could be resolved by carefully rescaling the system, but there is no systematic rule to determine a suitable rescaling.Second, working in multiple-precision offers a simple check on results: increase the precision, and confirm that the bounds U or L change very little.We have done this for all bounds reported here.Appendix A describes additional checks on numerical results, including how SDP computations could be part of fully rigorous computer-assisted proofs.
Example : deterministic upper bounds for the van der Pol limit cycle.To illustrate the application of Proposition 1, we have computed upper bounds ϕ = x 2 + y 2 ≤ U in the van der Pol system (1.3) for various 0.1 ≤ µ ≤ 5 and four different polynomial degrees of V u .Figure 2 shows the resulting U , along with estimates of ϕ obtained by integrating over the limit cycle using the fourth-order Runge-Kutta method with fixed time steps.As expected, increasing the degree of the storage function V u gives a better (smaller) upper bound U .For a given degree, the bounds become less tight as µ is raised because the shape of the limit cycle is more complicated (cf. Figure 1(b)).
The trivial lower bound 0 ≤ x 2 + y 2 is the best global lower bound possible for the van der Pol system.It cannot be improved because it is saturated by the trajectory staying at the unstable fixed point (x, y) = (0, 0).However, if one is interested only in the trajectories that tend to the stable limit cycle, then the lower bound of zero is not tight.Raising the lower bound requires a result that does not apply globally, but instead applies only to a subset of all possible trajectories.
Figure 2: Optimal upper bounds on ϕ = x 2 + y 2 for the van der Pol oscillator computed with the upper bound problem of (2.9) for different degrees of V u .The time average ϕ obtained by numerical integration (N.I.) of the system is also shown.

Bounds on local attractors
Often, one is interested in the average of a function ϕ(x) over only a subset of all possible trajectories, such as those tending to a particular attractor.If a bound on ϕ over such trajectories is to be tight, it should not apply to trajectories that tend to a different local attractor, nor should it apply to trajectories on unstable invariant structures that are not part of the attractor of interest (such as unstable fixed points, unstable limit cycles, or basin boundaries).However, global bounds of the kind derived using Proposition 1 are unlikely to be tight for a particular attractor since they must obey where x(t) is any trajectory of the dynamical system.Suppose we wish to bound possible values of ϕ, not for all trajectories but only for trajectories that eventually enter and remain inside a given absorbing domain T .To compute ϕ over any such trajectory, it suffices to begin averaging after the trajectory has permanently entered T , so dynamics outside of T can be ignored.This suggests the following modification of Proposition 1, where the inequality conditions are imposed only on T and, consequently, the resulting bounds are proven only for trajectories that permanently enter T .The comments made in Remark 1 still apply.
Proposition 2. Let ẋ = f (x) be a dynamical system with x ∈ R n , let ϕ(x) be a scalar function, let ϕ be its time average defined as in (1.1), and let T be an absorbing domain.If there exist differentiable functions V u (x), V l (x), and constants U , L such that then along any trajectory x(t) that permanently enters T , Applying Proposition 2 to a particular system requires specifying an appropriate absorbing domain.For instance, if one is interested in a local attractor A with basin B, it suffices to choose any T such that A ⊂ T ⊂ B. Bounds proven on T then apply to all trajectories that approach the attractor or are part of the attractor, but they need not apply to trajectories outside B.
Finding a mathematical description of T for a given attractor is difficult in general, but this too can be done using SoS programming [37,16].An example in which it is not difficult to choose T is when the trajectory to be excluded from the bound is a repelling fixed point.In this case it suffices to choose T = R n U for any small enough set U around that point.
If the absorbing domain can be specified as a semi-algebraic set-that is, defined by a set of polynomial inequalities and equalities-the conditions of Proposition 2 can be checked by SoS programming using the generalized S-procedure [36, Lemma 2.1].For instance, suppose T = {x | g(x) ≥ 0} for some polynomial g; proving a lower bound for trajectories entering T calls for SoS conditions that imply D l (x) ≥ 0 when g(x) ≥ 0. For this to be true, it suffices that there exist s(x) ≥ 0 such that D l (x) − s(x)g(x) ≥ 0. Strengthening these non-negativity constraints to SoS constraints and making similar arguments for the upper bound gives the following two SoS programs: where polynomial ansatze are specified for V u , V l , and s, and their free coefficients are the decision variables.Using similar ideas, the S-procedure can be generalized to semialgebraic T defined by multiple polynomial inequalities and equality constraints [36].
Example : deterministic lower bounds for the van der Pol limit cycle.Let us revisit the lower bound on x 2 + y 2 for the van der Pol oscillator.Suppose we want a bound that does not apply to the unstable fixed point at the origin but applies to the other trajectories, all of which approach the limit cycle.If such a lower bound is perfectly tight, it will equal the value of ϕ on the limit cycle.To apply Proposition 2, we must specify an absorbing domain T that contains the entire limit cycle but omits the origin.As described in Appendix B, the set T r = {(x, y) | g(x, y) = x 2 + y 2 − r 2 ≥ 0} is just such an absorbing domain for any r ≤ 1.
Figure 3 shows the lower bounds obtained by solving the lower bound program of (3.4) using the absorbing domains T 0.5 and T 1 .At a given polynomial degree, using the smaller Figure 3: Optimal lower bounds on ϕ = x 2 + y 2 over the limit cycle of the van der Pol system for different degrees of V l using the absorbing domains (a) The time average ϕ from numerical integration (N.I.) is also shown.absorbing domain T 1 gives a better bound.Results would likely be improved further by using T that closely approximate the attractor.This suggests a two-step procedure: using SoS techniques to find the smallest possible T around the attractor of interest [37,16], and then solving the SoS bounding programs (3.4).

Bounds for stochastic systems
Consider a stochastic dynamical system in which the random trajectories x(t) ∈ R n are bounded almost surely as t → ∞.The standard vector Wiener process ξ(t) ∈ R m is pre-multiplied by the matrix σ ∈ R n×m that describes the relative effect of each noise component ξ i on each state variable x j , and the overall noise strength is scaled by √ 2ε.We assume for simplicity that the noise is additive, meaning that σ does not vary with x, but the present analysis extends with only minor modifications to multiplicative noise in either the Itô or Stratonovich interpretation.We also assume that the system has reached statistical equilibrium, in which case the probability density of its trajectories, ρ(x) ≥ 0, decays at infinity and satisfies the stationary Fokker-Planck equation Suppose we wish to prove a constant lower bound ϕ ε ≥ L, where ϕ ε is the stationary expectation of ϕ(x) defined as in (1.2).We assume that such a stationary average exists, which is true for all ϕ(x) that don't grow too fast as |x| → ∞ (precise statements can be found in [27]).The subscript on ϕ ε indicates its dependence on the noise strength ε.
As in the deterministic case, our method of proof relies on a suitably chosen differentiable storage function V (x).For any V that does not grow too quickly as |x| → ∞, the stationary expectation ε∇ The third line above follows from the stationary Fokker-Planck equation (4.2).The second line follows from integration by parts, assuming that the boundary terms vanish-that is, assuming where ν(x) is the outwards unit normal to the sphere |x| = R, and dS is the surface element.The above condition holds in many cases, including any case where V is polynomial and ρ decays exponentially at infinity.When this condition holds, then so does the equality (4.4), hence the inequality (4.3) is equivalent to The above expectation cannot be evaluated without knowing ρ, but it is sufficient for the inequality to hold pointwise for all x.The lower bound ϕ ε ≥ L is therefore proven if we can find any differentiable V (x) satisfying the boundary integral condition (4.5) and the pointwise inequality The same argument with a reversed inequality sign gives a sufficient condition for an upper bound U on ϕ ε .We summarize these results in the following Proposition, which is the stochastic analog of Proposition 1.
be a stochastic dynamical system in a statistically stationary state with probability distribution ρ(x).If there exist differentiable functions V u , V l and constants U , L such that where D = σσ T , and if V u , V l grow slowly enough at infinity to each satisfy (4.5), then The inequality conditions (4.8a) and (4.8b) were derived by a different method in [8] for the case where D is the identity matrix (hence ∇ • D∇ = ∇ 2 ).The conditions differ from their deterministic counterparts (2.6a) and (2.6b), respectively, only in the addition of the diffusive terms ε∇ • (D∇V u ) and ε∇ • (D∇V l ).Similar results for non-negative ϕ were also stated in [13].
Like the deterministic bounds of § §2-3, the stochastic bounds and storage functions in Proposition 3 can be found numerically for a given ε if f and ϕ are polynomials.Letting V u and V l be polynomials also, we replace (4.8a) and (4.8b) with stronger SoS constraints to obtain the SoS programs min Since Proposition 3 relies on the boundary integral condition (4.5), the V u and V l constructed by the above SoS programs must satisfy this condition in order for U and L to be proven bounds.One way to guarantee this is to prove that ρ decays exponentially as |x| → ∞.
Example : stochastic bounds for the van der Pol oscillator.To illustrate the application of Proposition 3, we have bounded the stationary expectation ϕ ε = x 2 +y 2 ε for the stochastic van der Pol oscillator with µ = 1.For the σ matrix chosen above, the Wiener processes ξ acts on the y component alone.This corresponds to a stochastic physical force, as seen by rewriting (4.11) as a second-order equation for the position x(t), We have computed bounds using the SoS programs (4.10) for noise amplitudes 10 −6 ≤ ε ≤ 1 and three different polynomial degrees of V u , V l .Figure 4 shows the optimal bounds, along with several values of ϕ ε computed by numerically solving the stationary Fokker-Planck equation (4.2).Our solution of the Fokker-Planck equation employed finite differences with operator splitting as in [6], and the system was evolved to steady state by implicit Euler time stepping.The limiting expectation lim ε→0 ϕ ε equals the deterministic average ϕ on the limit cycle, reflecting the fact that the van der Pol system is stochastically stable.
The upper bounds U in Figure 4(a) are nearly tight at all noise strengths for V u of degree 10 and higher.They become less tight as ε increases, but they still correctly capture the increase in x 2 + y 2 ε that occurs when stochastic forcing "smears out" the van der Pol limit cycle.
The lower bounds L in Figure 4(b) are good for strong noise but not for weak noise.In fact, L → 0 as ε → 0 for any fixed polynomial degree of V l .To see that this is inevitable, observe that the diffusive term ε∇ • (D∇V l ) in (4.8b) vanishes as ε → 0 if ∇ • (D∇V l ) is bounded.As ε → 0, the stochastic SoS programs (4.10) reduce to the deterministic SoS programs (2.9), so the stochastic bounds can be no tighter than the global deterministic bounds.
In the van der Pol example, upper bounds are not affected by this phenomenon since the deterministic global upper bound on ϕ of §2 is sharp and ϕ ε → ϕ as ε → 0. However, a tight lower bound on ϕ ε at small ε is impeded by the unstable deterministic trajectory at the origin, which saturates the deterministic global lower bound ϕ ≥ 0. In the deterministic case, a tight local lower bound was achieved in §3 by ignoring part of phase space, but we cannot do so in the stochastic case since ρ > 0 everywhere.Instead, the diffusive term ε∇ • (D∇V l ) in the SoS constraint of (4.10) must remain O(1) near the origin as ε → 0. Clearly, this cannot be accomplished with polynomial V l .The next section illustrates how a non-polynomial V l can be used instead to improve the lower bounds on ϕ ε at small ε.

Bounds for stochastic systems with weak noise
An unstable trajectory of the deterministic system ẋ = f (x) can interfere not only with tight bounds on ϕ for a local attractor but also with tight bounds on ϕ ε in the corresponding stochastic system, at least when the noise is weak.This was illustrated in the previous section, where lower bounds on x 2 + y 2 ε in the van der Pol system were not tight: they approached zero as ε → 0 with polynomial V l of fixed degree.In this section we derive lower bounds on x 2 + y 2 ε that remain tight when ε is small.We do this using a non-polynomial V l that depends explicitly on ε.
The method described below applies to stochastic upper or lower bounds, provided that the deterministic trajectory interfering with the bounds is a repelling fixed point.For concreteness, suppose we are trying to prove a tight lower bound on an expectation ϕ ε that is strictly positive, as in the van der Pol example.Suppose also that the global lower bound on ϕ in the corresponding deterministic system, ẋ = f (x), is zero and is saturated by the repelling fixed point x = 0.For the lower bound on ϕ ε to remain larger than the deterministic lower bound as ε → 0, the stochastic bounding inequality (4.8b) must not reduce to its deterministic counterpart (2.6b).This requires that the diffusive term ε∇ • (D∇V l ) remains commensurate with the other terms in (4.8b), at least near the fixed point at the origin.
The term ε∇ • (D∇V l ) can remain O(1) near the origin as ε → 0 if V l develops a boundary layer there.Purely polynomial V l can have no such boundary layer, so we add a non-polynomial term and let ( Here, P (x) is a polynomial, α is a tunable constant, and ζ(x) is a positive definite quadratic form, for a suitable symmetric, positive definite matrix Z (denoted Z ≻ 0) to be determined.When ε is small, the logarithmic term dominates V l near x = 0, but P (x) remains significant away from the origin.Despite not being polynomial, the ansatz (5.1) is suitable to derive a SoS optimization problem for the bound because the relevant derivatives of V l are rational functions, Substituting these expressions into the bounding inequality (4.8b) gives a rational inequality with the positive denominator (ε + ζ) 2 .Multiplying this rational inequality by the denominator we obtain the equivalent polynomial inequality where ) Consequently, a lower bound L on ϕ ε can be calculated for a fixed, small ε by solving the optimization problem max (5.6) The above SoS program can produce lower bounds larger than zero for very weak noise strengths, meaning that these bounds are not spoiled by the fixed point at the origin.In fact, while the SoS bounding program (4.10) for polynomial V l reduces to its deterministic counterpart (2.9) as ε → 0, the above program does not.Instead, the polynomial L reduces to L 0 , retaining the term αζ (f • ∇ζ) that does not appear in the deterministic program.This term is due to the boundary layer built into V l around the origin.
Choosing the quadratic form ζ. If the coefficients of the positive quadratic form ζ (equivalently, the entries of Z) are tuned at the same time as the other coefficients, then the SoS constraint L ∈ Σ is quadratic in the decision variables, as opposed to linear.This produces a non-convex optimization problem for L that is harder than a standard SDP and requires a bilinear solver.Here we keep the problem convex by fixing a non-optimal ζ in advance and tuning only the other coefficients in V l .Even a non-optimal ζ provides a good lower bound at small ε provided that, for reasons explained in §6, we have not only ζ > 0 but also α ζ > 0 near the unstable fixed point.Quadratic ζ satisfying these two conditions can always be found if the point is a repeller but not if it is a saddle, which is why we have restricted ourselves to repellers.Further discussion of how to choose or optimize ζ is given in §6 and Appendix D.
Example : stochastic van der Pol oscillator with weak noise.We have used the methods of this section to derive lower bounds on x 2 + y 2 ε for the stochastic van der Pol oscillator (4.11).In the logarithmic ansatz (5.1) for V l , we used the positive quadratic form ζ = x 2 − xy + y 2 that satisfies ζ > 0 near the origin.Figure 5 shows the resulting  lower bounds computed using the SoS program (5.6) for various degrees of the polynomial P .When the noise amplitude is moderate or small, these bounds dramatically improve on the bounds of Figure 4(b) produced using purely polynomial V l .When ε 10 −3 , L is indistinguishable from the deterministic local lower bounds found in §3.

Bounds on local attractors: the vanishing noise limit
In this section we propose a second method for bounding deterministic averages on a local attractor.While in §3 we did this by segmenting the phase space, here we do it by adding noise.Many dynamical systems are stochastically stable, in the sense that the vanishing noise limit of a stationary expectation, lim ε→0 ϕ ε , is equal to the deterministic average ϕ on a particular attractor [42,9].This is true in the van der Pol example, where the vanishing noise limit of x 2 + y 2 ε is equal to x 2 + y 2 on the limit cycle.Such correspondence can be exploited to bound ϕ on a local attractor: by proving bounds on ϕ ε that hold as ε → 0, we obtain bounds that hold also for ϕ.In essence, the vanishing noise limit preserves the attractor while destroying unstable invariant structures.This idea was proposed in [8], though here we extend it by treating the ε → 0 limit rigorously, rather than numerically.
We continue the analysis of §5, still assuming that x = 0 is a repelling fixed point, ϕ(0) = 0, and ϕ > 0 on the attractor of interest.We suppose also that lim ε→0 ϕ ε = ϕ, meaning that the zero-noise limit of the expectation converges to the deterministic time average.Rigorous statements about when this property holds can be found in [42,9] and references therein.
Under these assumptions, we can obtain a lower bound L ≤ ϕ that is tight for the local attractor by deriving a lower bound L ≤ lim ε→0 ϕ ε .To this end, we recall that although the polynomial inequality L ≥ 0 considered in §5 suffices to prove the lower bound, we are really interested in proving the integral inequality (4.6).It is proven in Appendix C, under mild assumptions on the stationary distribution ρ, that (4.6) holds in the limit ε → 0 if there exists any γ > 0 such that L 0 ≥ γζ 2 .Dividing this relation by ζ (which is constructed to be positive definite) and rearranging gives the sufficient condition Since γ > 0 decreases L by an arbitrarily small amount, we can let the bound be implied by a strict inequality and set γ = 0.The lower bound we seek, L < lim ε→0 ϕ ε , is thus returned by the SoS optimization problem max where f = J 0 x denotes the linearized dynamics near the origin.This condition is needed because the SoS constraint becomes   near the origin.The only way the above condition can be satisfied for L larger than zero is if its first term is positive.Appendix D gives a way of constructing an admissible ζ when the unstable fixed point is repelling.When the unstable fixed point is a saddle, it is not generally possible to satisfy (6.3) for positive definite ζ.
Example : deterministic bounds for the van der Pol limit cycle-vanishing noise formulation.To demonstrate the methods of this section, we have computed lower bounds on the vanishing noise limit of x 2 + y 2 ε for the van der Pol system using the SoS program (6.2).The results serve also as deterministic lower bounds on ϕ for all trajectories approaching the limit cycle (strictly speaking, this conclusion requires proving that lim ε→0 x 2 + y 2 ǫ = x 2 + y 2 , which is outside our present scope).Bounds of the latter type appear also in §3, computed differently using the SoS formulation (3.4). Figure 6 shows the bounds obtained by fixing the degree of the polynomial P to 12 for the three different ζ defined in Table 1.The condition (6.3) is satisfied with α > 0 by ζ 2 and ζ 3 when µ = 2, and by The failure of the condition is why the bound using ζ 1 is poor for µ ≤ 4 − 2 √ 3 ≈ 0.54, while the bounds using ζ 2 or ζ 3 are poor near µ = 2. Taking the best of the three lower bounds at each µ gives a lower bound on x 2 + y 2 within 1% of the value found by numerical integration over the entire range 0.1 ≤ µ ≤ 5.

Future directions
Despite our success in computing bounds for the van der Pol oscillator, the study of more complicated dynamical systems presents several obstacles that require further investigations.We discuss some of these below, and put forward some preliminary ideas for future improvements.

Bounds for high-dimensional systems
The techniques we have developed in this paper can in principle be applied to systems of arbitrarily high but finite dimension.However, computing tight bounds typically requires the use of polynomial storage functions of large degree.Since the size of the SDP required to determine the existence of a SoS decomposition for a polynomial of degree d in n variables grows as n+d d [33,Theorem 3.3], relatively large SDPs must be solved even for a low-dimensional system such as the van der Pol oscillator.Moreover, the SDPs seem to be consistently ill-conditioned, and more so for larger systems.Our present methods are therefore practical only for systems of moderate dimension.
One way to improve the numerical conditioning of the SDPs is to try to rescale the system.A preliminary investigation showed that rescaling the van der Pol system such that the limit cycle is contained in the box [−1, 1] 2 dramatically improves the numerical conditioning, and bounds very similar to those reported above can be obtained without the need for multiple precision solvers.Researchers studying other aspects of dynamics using SDPs have had similar success rescaling the relevant dynamics to lie in [−1, 1] n .However, the appropriate rescaling is not generally clear a priori.
There are several possible ways to reduce the computational cost of constructing bounds.One option is to exploit special structures in the SDPs such as symmetries or sparsity (e.g.[11,28]).It might also be profitable to exploit special structure of the underlying dynamics, such as conserved quantities.Finally, it is sometimes possible to work with a lower-dimensional truncation of a high-dimensional system and bound the errors introduced by the truncation (e.g.[14,8,18]).

Bounds for systems governed by partial differential equations
The methods developed here apply only to finite-dimensional systems, and more work is needed to extend them to partial differential equations (PDEs).One idea is to project the PDE variables onto a finite Galerkin basis.Bounds can then be constructed using the finite-dimensional system, provided that the influence of the unprojected component can be controlled.This approach has been used successfully for nonlinear stability analysis of a fluid dynamical PDE [14,8,18].A second idea is to bound integrals of the PDE variables directly using the dissipation inequalities proposed in [1].In essence, these inequalities are integral constraints that replace the polynomial constraints of Proposition 1 when a finite-dimensional system is replaced by a PDE.Despite recent advances in the numerical implementation of integral inequalities [41], however, this technique has so far been applied only to simple examples.

Tight bounds for systems with saddle points
If deterministic bounds are spoiled by a saddle point x s that is not embedded in the attractor of interest, the S-procedure can be used to remove a set containing x s that is disjoint from the attractor.However, it can be difficult to establish that x s is indeed separate from the attractor.If x s is embedded in the attractor, alternative methods must be found.The vanishing-noise approach of § §5-6 seems promising, but the present formulation works only for repelling points, not saddle points.This is because (6.4) requires that the polynomial αζ(x) increases along all trajectories near the unstable fixed point, but αζ is monotonic and cannot increase along trajectories on the stable and unstable manifolds of x s simultaneously.
Similarly, the logarithmic ansatz (5.1) is not suitable when stochastic bounds with weak (but finite) noise are spoiled by a saddle point.In fact, the term f • ∇V in each of the inequalities in (4.8) has opposite signs along the stable and unstable manifolds of x s unless ∇V changes rapidly, and any negative contribution must be balanced by large Laplacians.This requires polynomials of large degree, making the SoS programs intractable in practice.
One possible solution is to find an ansazt for V that, similarly to (5.1) for repelling fixed points, lets stochastic bounds stay tight as ε → 0 by developing appropriate boundary layers.The following proposition suggests a way to deduce the ideal scaling of V in this limit.
Proposition 4. Consider the inequality that is a sufficient condition for the lower bound L ≤ ϕ ε .If S(x) is continuous and the stationary distribution ρ(x) is piecewise continuous, then the above inequality can be satisfied for the perfect lower bound L = ϕ ε only if S(x) ≡ 0 wherever ρ(x) > 0.
This statement is proven by adding ϕ − L ε = 0, which holds because L is a perfect bound, and expression (4.4) to find S ε = 0. Since S(x) is continuous and non-negative, S ε = 0 is possible only if S(x) = 0 wherever ρ(x) > 0, as claimed.An analogous statement holds for the perfect upper bound.
Proposition 4 suggests that as V l is improved and L gets closer to its perfect value of ϕ ε , the inequality S(x) ≥ 0 gets closer in some sense to being an equality.This motivates us to consider solutions V (x) to the equation where L = ϕ ε and ρ obeys the decay condition (4.5).A solution exists only when L = ϕ ε , as seen by taking the expectation of the equality and using (4.4).Asymptotic analysis of V (x) in the above equality as ε → 0 might suggests a good ansatz for V l in the corresponding inequality.

Assessing the quality of numerical bounds
Proposition 4 also provides a good way of assessing a posteriori whether the stochastic bounds obtained with SoS optimisation are nearly sharp.This is useful, for instance, when data from numerical simulations are not available for comparison.For example, if a lower bound is close to the exact value of ϕ ε , V l should be close to the solution of (7.2) whenever ρ(x) = 0.In addition, an argument similar to the proof of Proposition 4 for deterministic bounds shows that if an upper or lower bound is perfect, then for any trajectory x(t) on the attractor.If a good lower bound is proven using V l , for instance, then Vl should be very close to ϕ − ϕ(x) on the attractor.
For illustration, we have solved the lower bound program (6.2) for the van der Pol oscillator with P (x) of degree 6 and degree 12.The degree-6 lower bound of L = 2.74 is rather poor, while the degree-12 bound of L = 4.11 is nearly tight.Figure 7 shows the variation of Vl (t) along the limit cycle in each case and compares it to the value of ϕ−ϕ(x) along the limit cycle.The Vl (t) and ϕ − ϕ(x) curves match very closely in the degree-12 case but not in the degree-6 case.).The logarithmic ansatz (5.1) for V l was used with P of order 6 and 12.The storage functions were calculated with (6.2), while the limit cycle is from numerical integration.

Formulating computer-assisted proofs
The bounds we have computed for the van der Pol system all rely on the numerical solution of SDPs, so they are subject to roundoff error.It is unlikely that roundoff error has led to any invalid bounds in our calculations because we have solved the SDPs using multipleprecision arithmetic and checked the results by increasing the precision.Nonetheless, such computations cannot be considered rigorous to the standard of computer-assisted proofs because of the presence of roundoff errors.Two different methods have been proposed for obtaining rigorous results from numerical SDP solutions.The first method is to project an approximate numerical solution onto an exact solution in terms of rational numbers [34,20].The second method is to use perturbation analysis, made rigorous by interval arithmetic, to construct a small interval around the approximate numerical solution in which the exact solution is guaranteed to lie.This approach has been implemented in the software VSDP [15], but the relevant SDPs would need to be constructed manually since there is no sum-of-squares parser available that incorporates interval arithmetic.

Conclusions
In this work we have presented computer-assisted methods for deriving bounds on average quantities in both deterministic and stochastic dynamical systems using sum-of-squares programming.We have given particular attention to proving bounds that apply only to trajectories approaching a particular attractor.One method is to use the S-procedure to omit segments of phase space on which the bounds do not need to hold.Another strategy is to remove unstable invariant structures by adding noise to the system.This idea, proposed previously in [8], has been extended here by analyzing the weak and vanishing noise cases when the unstable structures to be omitted are repelling fixed points.These methods give improved bounds for weak but finite noise, and also in the rigorous limit of vanishingly weak noise.Our methods have worked well when applied to the van der Pol oscillator.The best deterministic and stochastic bounds proven throughout are summarized in Figure 8.In the deterministic case, we obtained upper and lower bounds on the infinite-time average of x 2 + y 2 over the limit cycle that are all within 1% of the "true" values found by numerical integration.In the stochastic case, bounding the stationary expectation of x 2 + y 2 at a  variety of noise strengths, we obtained upper bounds all within 1% and lower bounds all within 10% of the "true" values found by solving the Fokker-Planck equation numerically.
Moving to dynamical systems other than the van der Pol oscillator, whether the methods we have described can yield similarly tight deterministic and stochastic bounds depends on the details of those systems.If the dimension of a system is small enough for SoS optimization to be computationally feasible, we expect that tight global bounds on deterministic averages can be obtained using the methods of §2.In fact, since the first draft of this work, the bounding techniques we have presented for deterministic systems have been successfully applied to the design of control systems for fluid flows [19,23].With stochastic forcing that is not too weak, we expect that fairly tight bounds on stationary expectations also can be obtained using the methods of §4.As the noise strength decreases, the tendency of unstable invariant structures to spoil tight bounds can be combatted by the methods of §5 if these structures are repelling fixed points.New methods will be needed to maintain tight bounds as noise becomes weak in systems with other unstable structures, including saddle points.Lastly, the methods of §3 can be used to bound averages on trajectories within an absorbing domain and, in particular, on a single attractor.This much is true for any dynamical system, but whether these bounds can be tight depends on the nature of the attractor.For instance, bounds for a chaotic attractor must apply not just to generic chaotic trajectories but also to all trajectories embedded within the attractor, such as saddle points and unstable orbits.What can be deduced about chaotic attractors is of particular interest if the methods presented here are to be applied to complex systems of physical and engineering relevance.

A Towards rigorous bounds
As explained in §7.5, techniques are available to control the roundoff errors in the solution of the SDPs and compute rigorous bounds.However, these may be impractical for large problems.We are also unaware of any parsers for SoS programs that incorporate rigorous computations.For these reasons, we illustrate a method to produce rigorous bounds that can be implemented with existing software.For definiteness, we consider the upper bound problem in (2.9).Comparison between the optimal deterministic upper bounds on ϕ = x 2 + y 2 of §2 and the certified bounds obtained with Algorithm 1 for V u of degree 10.The "true" value of ϕ from numerical integration (N.I.) is also shown.

B Construction of an absorbing domain for the van der Pol oscillator
To show that T r = {(x, y) | x 2 + y 2 − r 2 ≥ 0} is an absorbing domain for the van der Pol oscillator for r ≤ 1, let us reverse the direction of time in (1.3) and consider the energy E = x 2 + y 2 of the system ẋ = −f (x), for which the origin is a stable fixed point.meaning that E ≤ 0 when |x| ≤ 1.Therefore, any contour of E which is contained in the strip |x| ≤ 1 defines the boundary of a region of attraction of the origin for the timereversed oscillator.One concludes that in the original system all orbits (except of course the fixed point at the origin) will leave the ball x 2 + y 2 < r 2 if r ≤ 1, i.e.T r is an absorbing domain.

C Bounds with vanishing noise
To prove a lower bound on ϕ ε we need to prove inequality (4.6).The ansatz (5.1) for V allows us to rewrite this condition as where L is as in (5.4).Since ζ is quadratic in x, it can be verified that L 0 is the dominant term in L for any fixed x = 0 in the limit ε → 0.Moreover, L 0 is the dominant term as |x| → ∞ because it contains the monomials of highest order.The condition L 0 ≥ γζ 2 then implies that L(x) ≥ 0 in the limit ε → 0 for all x, with the possible exception of a ball B R of radius R ∼ ε η , 0 < η < 1/2, where the integrand in (C.1) becomes singular.
For definiteness, we will write R = Cε η for some constant C. Thus, we conclude that Let us proceed term by term and let us assume that ρ is bounded in B R uniformly as ε → 0. This is a reasonable assumption, since we are assuming that the point x = 0 is a repeller and hence does not belong to the attractor of the deterministic system on which we require the bound.Since P and ϕ are continuous,

D Construction of ζ
If the unstable fixed point is a repeller and J 0 can be diagonalized, a quadratic form ζ that satisfies (6.4) can be constructed from the eigenvectors of J 0 .Let U denote the matrix of eigenvectors of J 0 and Λ be the usual diagonal matrix of eigenvalues, such that Note that Λ is a positive definite matrix since the unstable point is repelling.Letting w = U −1 x, an appropriate choice of ζ is In fact, we have f which is positive for x = 0 because Λ is positive definite.We conclude that (6.4) holds for α > 0. Note, however, that this is not the only possible choice of ζ.

Figure 1 :
Figure 1: (a) Phase portrait of the van der Pol oscillator for µ = 1.The limit cycle and the unstable fixed point at (x, y) = (0, 0) are highlighted.(b) Limit cycles for various values of µ.

Figure 4 :
Figure 4: (a) Upper bounds and (b) lower bounds on ϕ ε = x 2 + y 2 ε for the stochastic van der Pol oscillator for µ = 1 as a function of the noise strength ε.The deterministic (ε = 0) average ϕ = x 2 + y 2 ≈ 4.118 and the expectation ϕ ε computed from the Fokker-Planck equation for selected values of ε are also shown.Improved lower bounds appear in Figure 5.

(6. 2 )
Like the finite-noise program (5.6), the above program can be made convex by fixing in advance the positive definite quadratic form ζ. The chosen ζ must satisfy α ζ > 0 near the origin.That is, it must satisfy α f • ∇ζ > 0,(6.3)

Figure 6 :
Figure6: Lower bounds on ϕ = x 2 + y 2 for the van der Pol oscillator obtained with the vanishing noise program (6.2) for P of degree 12, along with the values of ϕ found by numerical integration (N.I.).

Figure 7 :
Figure7: Comparison of Vl along the limit cycle of the deterministic van der Pol oscillator for µ = 1 (where ϕ = x 2 + y 2 ≈ 4.118).The logarithmic ansatz (5.1) for V l was used with P of order 6 and 12.The storage functions were calculated with (6.2), while the limit cycle is from numerical integration.

Figure 8 :
Figure 8: Best bounds obtained on averages of ϕ = x 2 + y 2 for the van der Pol oscillator.(a) Bounds on the deterministic time average ϕ over the limit cycle for various damping strengths µ.(b) Bounds on the stochastic expectation ϕ ε with µ = 1 and various noise strengths ε.

Figure 9 :
Figure 9: Comparison between the optimal deterministic upper bounds on ϕ = x 2 + y 2 of §2 and the certified bounds obtained with Algorithm 1 for V u of degree 10.The "true" value of ϕ from numerical integration (N.I.) is also shown.