## Abstract

The first part of this paper reviews the application of the sum-of-squares-of-polynomials technique to the problem of global stability of fluid flows. It describes the known approaches and the latest results, in particular, obtaining for a version of the rotating Couette flow a better stability range than the range given by the classic energy stability method. The second part of this paper describes new results and ideas, including a new method of obtaining bounds for time-averaged flow parameters illustrated with a model problem and a method of obtaining approximate bounds that are insensitive to unstable steady states and periodic orbits. It is proposed to use the bound on the energy dissipation rate as the cost functional in the design of flow control aimed at reducing turbulent drag.

## 1. Introduction

SOS stands for sum of squares. In recent years, it has become common for references to SOS of polynomials. It is used to refer to a recent discovery that the SOS decomposition for a polynomial can be computed via semidefinite programming (SDP) [1,2]. This technique provided a constructive method for generating Lyapunov functions for systems whose dynamics can be described by polynomial functions, and made many dynamical-system problems computationally tractable. This paper gives a review of the progress made so far in using this breakthrough in fluid dynamics, and continues by offering new ideas of using SOS in the areas of flow stability, bounds for time-averaged characteristics of turbulent flows and flow control.

### (a) The idea of sum of squares

The SOS technique is useful in situations in which one must verify that a given polynomial function of several variables is positive-definite, or when a positive-definite polynomial must be constructed subject to some algebraic constraints. The classical algebraic–geometry problem of verifying positive-definiteness of a general multi-variate polynomial is NP-hard. This means that, in general, its numerical solution is far beyond the capabilities of modern computers. SOS methods are based on the recognition that a sufficient condition for a polynomial function to be positive-definite is that it can be rewritten as an SOS of lower-degree polynomial functions. While the full technical details are complicated, the basic idea is simple [1,2]. An even-degree polynomial *P*(**a**) of variables (*a*_{1},…,*a*_{n}) can always be represented as
1.1
for some (non-unique) matrix *Q*, where **m** is a vector of monomials of **a**, that is, **m**=(*a*_{1},…,*a*_{n},*a*_{1}*a*_{1},…,*a*_{k}*a*_{l},…,*a*_{n}*a*_{n},…). Taking into account all the monomials of degree up to half of the degree of *P* is sufficient. Without loss of generality, the matrix *Q* can be assumed to be symmetric. If there exists some positive-semidefinite *Q* (denoted *Q*≽0) satisfying (1.1), then a linear transform of the variables can reduce it to a diagonal matrix with all the elements on the main diagonal being non-negative. Because a linear combination of monomials is a polynomial, this gives a representation of the polynomial *P* as an SOS of other polynomials.

The resulting problem is now a problem of finding a matrix *Q* subject to a set of linear constraints (1.1) and a constraint *Q*≽0. This is a feasibility problem on a set of convex constraints. One can also specify that *Q* should maximize a given linear function of *Q*. This is not needed for a Lyapunov stability analysis, but will be useful in other applications discussed later. The problem turns out to be reformulated as an optimization problem in the form of SDP. SDP optimization problems are convex and tractable (e.g. solvable in a number of operations that is a polynomial function of the problem size [3]), and theoretical and algorithmic research into solving such problems is extremely active [4,5]. A variety of well-supported software codes for solving such problems are freely available, both for SOS problems in particular [6] and for SDP problems in general [7,8]. More recently, research has focused on the efficient solution of SDP problems arising from SOS problems, with particular emphasis on robust optimization [9,10], exploitation of structure and sparsity [11–13] and decomposition techniques [14] in large problems. As a result of these advances, the SOS approach has found extensive applications in stability analysis, control theory and many other fields, including applications in aeronautics [15,16]. The use of SDP also has a long history in systems analysis and control, because many control problems can be formulated and solved this way [17].

The SOS approach has an immediate application in the stability theory of systems with polynomial dynamics, that is, the systems governed by ordinary differential equations with polynomial right-hand side **f**(**a**),
We assume that **f**(0)=0. A steady solution of such a system is stable if one can find a continuously differentiable Lyapunov function *V* (**a**) satisfying each of the following conditions [18], theorem 4.1:
1.2
1.3
1.4
There is no known method to construct such a function for an arbitrary system of nonlinear ordinary differential equations. However, if **f** is polynomial as in our case, and a polynomial Lyapunov function *V* is of interest, then (1.3) and (1.4) can be relaxed to SOS constraints, that is, requirements that *V* (**a**) and −∇_{a}*V* (**a**)⋅**f**(**a**) are SOS of polynomials. To accommodate the strict inequalities in (1.3) and (1.4), these conditions can be slightly modified by buffer functions (e.g. by replacing (1.3) with *V* (**a**)≥*ϵ*|**a**|^{2} with a very small *ϵ*). The SOS constraints are convex and reducible to the condition of matrices being positive-definite, as explained by (1.1), whereas (1.2) is linear. This makes the problem numerically tractable. For systems of relatively modest dimension, a polynomial Lyapunov function, if it exists, can be found numerically using the freely available Matlab toolboxes SOSTOOLS [19,20] or YALMIP [21,22], along with an SDP solver such as SEDUMI [7] or SDPT3 [8]. For the systems obtained by Galerkin approximation of the Navier–Stokes equations a calculation using about 10 modes does not require a supercomputer. The volume of the calculations grows fast with the system size, and special methods discussed in the following sections might be needed.

To illustrate how SOS can be used to analyse stability of a nonlinear system, we provide an example. Please refer to appendix A for more details.

Consider the system For this system, a quadratic Lyapunov function does not exist, as can be checked by considering a quadratic polynomial with unknown coefficients and calculating its time derivative. The same conclusion can be algorithmically tested using SOSTOOLS. Therefore, a simple ‘energy’ Lyapunov function structure in this case is inadequate. At the same time, global asymptotic stability of the zero equilibrium is difficult to demonstrate, as using a quartic Lyapunov function by hand becomes too involved.

Using the `findlyap` command in SOSTOOLS, we can search for a quartic Lyapunov function. The command returns a Lyapunov function of the form
The conditions (1.2)–(1.4) can be tested using the `findsos` command in SOSTOOLS. Demo 2 in the SOSTOOLS manual gives a worked example of how this toolbox can be used to construct Lyapunov functions in more complicated cases.

## 2. Applying sum of squares for stability analysis of fluid flows

Flow stability theory is a well-established part of fluid dynamics, with Helmholtz, Kelvin, Rayleigh and Reynolds among its pioneers [23]. The bulk of the work on hydrodynamic stability is concerned with infinitesimal perturbations. This allows one to represent the solution of the governing equations as a sum of the steady solution and a small perturbation, and to neglect the nonlinear terms. The resulting linear problem is much easier to solve. Linear stability of canonical flows is now largely a closed area of research, with the centre of gravity shifting to development of efficient numerical methods applicable to complex flows encountered in practice. Importantly, linear stability analysis can reveal instability but cannot prove stability. This is because steady flows are often stable with respect to infinitesimally small perturbations, but unstable with respect to finite perturbations. Moreover, the finite amplitude required to destabilize the flow is often small. Hence, it is the stability with respect to finite disturbances and the control of finite disturbances that represent the major practical interest.

Typically, flows are stable if the Reynolds number *Re* is small enough, but there is a value *Re*_{l} such that, if *Re*>*Re*_{l}, the flow is unstable with respect to arbitrary small perturbations. On the other hand, there is a critical Reynolds number *Re*_{c} such that, if *Re*<*Re*_{c}, the flow is stable with respect to disturbances of any amplitude. Finding *Re*_{c} is difficult because this is a nonlinear problem. Serrin [24] demonstrated that, for each incompressible flow in a closed domain, there is *Re*_{e} such that the energy of an arbitrary perturbation decreases monotonously if *Re*<*Re*_{e}. This, of course, means that the flow is globally stable. Moreover, Serrin demonstrated that the problem of determining *Re*_{e} can be reduced to a linear eigenvalue problem similar to the eigenvalue problems arising in linear stability theory. Naturally, *Re*_{e}≤*Re*_{c}≤*Re*_{l}. In many cases, the difference between *Re*_{e} and *Re*_{l} is large: for example, for a pipe flow . Therefore, while finding *Re*_{e} and *Re*_{l} is relatively easy, more powerful methods are required for estimating *Re*_{c}. Progress in understanding instability with respect to finite disturbances was made in the late 1980s in the works on non-modal stability and their extensions to nonlinear settings (see reviews [25,26]), but no methods for determining *Re*_{c} have yet emerged from these works. SOS is the natural way forward in the analysis of stability with respect to perturbations of arbitrary amplitude.

### (a) Direct application of the sum of squares technique to Navier–Stokes equations

Fluid flows are governed by partial differential equations. In the case of incompressible flow, the Navier–Stokes equations have the form
2.1
with appropriate boundary conditions. A steady solution , of the system (2.1) is globally stable if, for each *ϵ*>0, there exists some *δ*>0 such that at time *t*_{0} implies that for all time *t*≥*t*_{0}. It is globally asymptotically stable if in addition as time for any initial conditions. In analogy with a Lyapunov function, one can introduce a Lyapunov functional *V* [**u**′]>0 ∀**u**≠0, *V* [0]=0, , where . With **u** being a particular solution of (2.1), *V* becomes a function of time. If d*V*/d*t*<0 for any solution except , then the steady solution is globally stable.

To illustrate the main difficulties in applying the SOS technique to (2.1) and the ideas of overcoming them, consider the candidate Lyapunov functional in the form of a volume integral
where *P* is a polynomial and *X* is the flow domain. (Polynomial *P* can also depend on the derivatives of **u**′, which complicates the matter and is an area of intensive research.) Then
One can then substitute ∂**u**′/∂*t* from (2.1) to obtain
where *K* is a polynomial function of **u**′,∇**u**′,∇^{2}**u**′ and ∇*p*, in which the physical coordinate vector **x** is also involved via **F** and . The full expression for the kernel *K* is too involved to be given here. Now, if it were possible to find *P*(**u**′) such that *V* is positive-definite and *K*<0 everywhere in *X* for any values of its arguments, then *V* would be proved to be a Lyapunov functional. However, this is hardly possible. The idea [27], in general applicable to any partial differential equations and not only to (2.1), is to seek generic identities of the form
2.2
where polynomials *A*_{i} contain arbitrary coefficients. Typically, such identities can be obtained by integration by parts. For example, for a flow in a closed domain with a no-slip condition on the wall
for any scalar *α*(**u**′,∇**u**′,∇^{2}**u**′,∇*p*,**x**) and any solenoidal **u**′, provided that the normal component of **u**′*α* at the boundary of *X* is zero. Then, the requirement *K*<0 can be relaxed to , and the SOS approach can be tried to tune the free coefficients in *P* and *A*_{i} to make a sum of squares. This strategy proved to be successful for some simple partial differential equations. It appears that there is only one reported attempt to apply this approach to the Navier–Stokes equations [28,29]. However, in [29], the results were reported only for *P*=∥**u**′∥^{2}, so that the Lyapunov functional was equivalent to the Lyapunov functional in the standard energy stability approach, but, even in this case, it appears that the highest *Re* for which the stability was proved was well below the known [30] energy stability limit for the same problem of Hagen–Poiseuille flow with two-dimensional perturbations. In [28], for the same problem, results with a different *P* were also reported. However, in this generic case, *K* contains a pressure term, and the procedure described in [28] does not explain how this term was treated; moreover, the results for *P*≠∥**u**′∥^{2} were not mentioned in the later paper [29] of the same authors. In any case, the identities of type (2.2) found in [28,29] for the Hagen–Poiseuille flow are of definite interest.

Unlike the SOS part of the approach, which can be done by a computer, there is no systematic procedure for finding identities (2.2) for a particular flow. This, therefore, remains a challenge for the ingenuity of the researcher. We now describe an alternative, systematic, method of applying SOS to the Navier–Stokes equations, which does not require identities (2.2).

### (b) Reducing the Navier–Stokes equations to an uncertain finite-dimensional dynamical system

Theoretically, a straightforward way of applying the SOS approach to the Navier–Stokes equations would be to use a truncated Galerkin expansion in order to approximate the full Navier–Stokes equations with a finite-dimensional system of ordinary differential equations. In practice, in order to achieve good approximation, the dimension of the system must be sufficiently large. Then, the resulting SOS problem turns out to be too large for contemporary SDP algorithms running on contemporary computers.

The idea for overcoming this difficulty was proposed in [31,32]. The perturbation velocity is partitioned into a sum of a projection onto a finite-dimensional subspace and a residual infinite-dimensional velocity field,
2.3
where the finite Galerkin basis **u**_{i}, *i*=1,…,*k* is an orthonormal set of solenoidal vector fields with appropriate inner product, the residual perturbation velocity **u**_{s} is solenoidal and orthogonal to all the **u**_{i}, and both **u**_{i} and **u**_{s} satisfy the boundary condition imposed on **u**.

Substituting (2.3) into the Navier–Stokes equations and projecting them onto the basis **u**_{i} gives
2.4
where **a**=(*a*_{1},…,*a*_{k})^{T}. Explicit expressions for **f** and ** Θ**=

*Θ*_{a}+

*Θ*_{b}+

*Θ*_{c}in terms of

**a**,

**u**

_{s},

**u**

_{i}and are [32]: ,

*Θ*

_{ai}(

**u**

_{s})=〈

**u**

_{s},

**g**

_{i}〉,

*Θ*

_{bi}(

**u**

_{s},

**a**)=〈

**u**

_{s},

**h**

_{ij}〉

*a*

_{j},

*Θ*

_{ci}(

**u**

_{s})=〈

**u**

_{s},

**u**

_{s}⋅∇

**u**

_{i}〉,

*L*

_{ij}=(1/

*Re*)〈

**u**

_{i},∇

^{2}

**u**

_{j}〉+〈

**u**

_{i},

*A*

**u**

_{j}〉, The inner product 〈

**w**

_{1},

**w**

_{2}〉 is the integral of

**w**

_{1}⋅

**w**

_{2}over the flow domain, and

**h**

_{ij}=

**u**

_{j}⋅∇

**u**

_{i}−

**u**

_{i}⋅∇

^{T}

**u**

_{j}.

Note that ** Θ**=0 for

**u**

_{s}=0, and for

**u**

_{s}=0 equation (2.4) is a usual truncated Galerkin approximation of the Navier–Stokes equations written in velocity perturbation, of the type commonly used in numerical calculations. With

**u**

_{s}≠0, however, (2.4) is not a closed system of equations.

Instead of considering in full the dynamics of **u**_{s}, which is described by a system of partial differential equations, it was proposed [32] to find a bound for ** Θ**. The energy of the residual velocity field

*q*

^{2}=∥

**u**

_{s}∥

^{2}/2 satisfies the equation 2.5 where explicit expressions for

*Γ*(

**u**

_{s}) and

*χ*(

**u**

_{s},

**a**) are also available [32]. If

**u**

_{i}are chosen to be the eigenfunctions of the energy stability problem, then

*χ*(

**u**

_{s},

**a**)=0 identically. We assume that

**u**

_{i}are the eigenfunctions of the energy stability problem. Then where

*κ*

_{s}is the largest eigenvalue of all the energy stability problem eigenvalues excluding those that correspond to

**u**

_{i}. By properly choosing

*k*and

**u**

_{i}, it is always possible to ensure that

*κ*

_{s}<0. A bound in the form |

**|**

*Θ*^{2}≤

*p*(

**a**,

*q*

^{2})=

*c*

_{1}|

**a**|

^{2}

*q*

^{2}+

*c*

_{2}

*q*

^{4}, where |⋅| represents the standard Euclidean norm in , is also available [32], and an important improvement described in [33] gives a simple method of calculating the constants

*c*

_{1}and

*c*

_{2}.

The next step is to abandon **u**_{s} entirely, keeping only the bounds in terms of **a** and *q*^{2} on the functions of **u**_{s}. This results in the uncertain dynamical system
2.6a
2.6b
2.6c
This finite-dimensional system is uncertain in the sense that ** Θ** and

*Γ*can vary arbitrarily within the constraints imposed by the inequalities. Accordingly, its solution is not unique. However, all possible solutions of the original Navier–Stokes system are within the solutions of (2.6), and therefore any statements valid for arbitrary solutions of (2.6) hold true for the solutions of the Navier–Stokes equations.

### (c) Sum of squares analysis of stability of fluid flows

#### (i) Model example for a truncated system

The theoretical background of the approach outlined in §2*b* was given in [32] in full, including the derivation of the uncertain system and the bounds, whereas the numerical example was given for the truncated system only, that is, for (2.6a) with ** Θ**=0. Nevertheless, this example is of interest, because it highlighted a certain difficulty in the stability analysis of fluid dynamical systems and a way of overcoming it, and showed by how much the global stability bound might be improved through the SOS approach. For the truncated system, the Lyapunov condition requires that the Lyapunov function is positive-definite,

*V*(

**a**)≻0, and that its time derivative along the system trajectory is negative-definite, . For using the SOS approach,

*V*(

**a**) has to be sought in the form of a polynomial, and a positive-definite polynomial has to be of an even degree. Note now that because the nonlinearity of the incompressible Navier–Stokes equations is quadratic,

**f**is a quadratic polynomial. Then, with

*V*of even degree, in general case

**f**⋅∇

_{a}

*V*is of an odd degree and cannot, therefore, be sign-definite. This difficulty can be overcome if one takes into account that the nonlinear terms of the Navier–Stokes equations are energy-conserving. For the Galerkin basis orthogonal in

*L*

_{2}, as is the basis we consider, it means that

**f**⋅∇

_{a}

*E*, where

*E*=|

**a**|

^{2}, is a polynomial not of the third but only of the second degree. This allows one to search for a Lyapunov function of the form

*V*=

*E*

^{n}+

*p*

_{2n−1}(

**a**), where

*p*

_{2n−1}is a polynomial of degree 2

*n*−1.

The SOS approach was applied to the ninth-order model [34] of the Couette flow, that is, a shear flow of a fluid between two infinite parallel plates. This model was specifically designed to reproduce all the major features of the flow, including the transition to turbulence, streaks and vortices, periodic orbits, chaotic sets, and so on [35]. For this system, the energy stability approach gives *Re*_{e}=7.5. Using the SOS approach, the largest Reynolds number for which a Lyapunov function was found was *Re*_{SOS}=54.1 [32]. Numerical results [35] suggest that the system becomes unstable at *Re*≈80. This implies that the SOS approach can be significantly more effective than the energy stability approach.

#### (ii) Sum of squares analysis of the stability of an uncertain system

To prove the global stability of the uncertain system (2.6) using SOS, one can attempt to find a positive-definite polynomial Lyapunov function *V* =*V* (**a**,*q*^{2})≻0 such that
2.7
for all **a**≠0 and *Γ* and ** Θ** satisfying (2.6c). If one chooses as a candidate Lyapunov function
2.8
then
and the situation reduces to the usual global stability condition using energy functions. Note that this is the case when energy stability eigenfunctions are used as the Galerkin basis in the SOS analysis, which is the assumption we made earlier. Consequently, if energy can be used as a Lyapunov functional for the Navier–Stokes equations for some Reynolds number

*Re*, then the choice (2.8) will satisfy the conditions (2.7).

When the system remains globally stable for Reynolds numbers beyond the energy stability limit, there exists a polynomial in (*a*,*q*^{2}) that serves as a Lyapunov function at least for *Re* just beyond *Re*_{e}. A constructive proof is given in [32]. However, not every positive-definite polynomial can be represented as a sum of squares.

It remains to convert (2.7) to an SOS problem. Assuming for the sake of simplicity that ∂*V*/∂*q*^{2}>0 and using (2.6c), (2.7) can be rewritten as
2.9
Note that this involves ** Θ** while the polynomial bound (2.6c) is a bound for |

**|**

*Θ*^{2}. Hence, an additional step is required. Using the Schwarz inequality and (2.6c) gives a sufficient condition for satisfaction of the inequalities (2.7) and (2.9), 2.10 One can now introduce an additional scalar variable

*z*

_{0}and a vector variable

**z**with components

*z*

_{1},…,

*z*

_{n}and define a polynomial It is relatively straightforward to verify that (2.10) is equivalent to 2.11 The advantage of (2.11) over (2.10) is that

*W*(

*z*

_{0},

**z**,

**a**,

*q*

^{2}) is a polynomial, whereas (2.10) involves a square root of the polynomial

*p*(

**a**,

*q*

^{2}). This advantage was achieved, however, at the expense of doubling the number of independent variables. In [32], an equivalent reduction was performed in the general case when the Galerkin basis does not necessarily consist of the eigenfunctions of the energy stability problem.

In summary, to prove the global stability of the uncertain system (2.6), and, hence, the global stability of the flow, for which the uncertain system was derived, it is sufficient to find a positive-definite polynomial Lyapunov function *V* (**a**,*q*^{2}) such that both ∂*V*/∂*q*^{2} and *W*(*z*_{0},**z**,**a**,*q*^{2}) are also positive-definite. Importantly, the coefficients of the polynomials ∂*V*/∂*q*^{2} and *W*(*z*_{0},**z**,**a**,*q*^{2}) are linear functions of the coefficients of *V* . This allows one to compute constraint-admissible coefficients of *V* , if they exist, using the existing packages SOSTOOLS [6] or YALMIP [21,22].

#### (iii) Application to a version of rotating Couette flow

The first application of this approach to a fluid flow is described in [33,36]. The rotating Couette flow is a well-known benchmark for flow stability studies [37]. It is a flow between two rotating co-axial cylinders. Assuming that the gap between the cylinders is much smaller than the cylinder radius, a Cartesian coordinate system **x**:=(*x*,*y*,*z*) can be introduced in such a way that the axis of rotation is parallel to the *z*-axis, whereas the circumferential direction corresponds to the *x*-axis. Only flows independent of *x* were considered in [33,36]. The flow velocity was represented as (*y*+*u*′,*v*′,*w*′), so that **u**′=(*u*′,*v*′,*w*′) was the velocity perturbation. Under these assumptions, the governing equations for **u**′ are [37]
where *Ω* is a non-dimensional parameter characterizing the Coriolis force, and *p*′ is pressure. For the sake of simplicity, the flow is assumed to be independent of *x* and 2*π*-periodic in *y* and *z*, *u*′ and *v*′ are assumed to be odd in *y* and even in *z*, whereas *w*′ is odd in *z* and even in *y*.

This system is linearly unstable when , 0<*Ω*<1. Perturbation energy decays monotonously if . Note that similarly to the case of the classic rotating Couette flow formulation, the energy stability limit is independent of *Ω*, and that the energy stability limit and the linear stability limit coincide for *Ω*=1/2. It might reasonably be expected [36] that, in contrast to the classic rotating Couette flow formulation, this system is globally stable if it is linearly stable owing to the symmetry constraints imposed on it. For a particular selection of six energy stability eigenfunctions used as the Galerkin basis it is shown in [36] that, for *Ω*∈(0.253,0.747), the resulting uncertain system is globally stable if the flow is linearly stable. For *Ω*∈(0,0.253]∪[0.747,1), the uncertain system was shown to be globally stable for , and globally unstable for larger *Re*. Because the absence of global stability of the uncertain system does not imply that the true flow is globally unstable, one can expect that the result can be improved by increasing the number of modes explicitly accounted for in the uncertain system. However, an attempt to do this was not successful [36]. This might imply that the quality of the uncertainty bounds derived and used in [32,33,36] deteriorates as the dimension of the uncertain system increases.

Most importantly, the results showed that the range of *Re* in which the flow has been proved to be globally stable goes beyond the energy stability limit.

## 3. A look ahead

Here, several new ideas are introduced, and some indicative new results are presented, revolving around the application of the SOS approach to fluid dynamics problems other than flow stability.

### (a) Energy dissipation bounds

While obtaining rigorous solutions of the full Navier–Stokes equations for the turbulent flow regime appears to be impossible, rigorous bounds on various characteristics of turbulent flows can be derived. Bounds on time-averaged momentum transport were obtained in [38,39]. Bounds on the time-averaged energy dissipation rate were obtained in [40], using the background-flow idea. These studies were based on optimizing quadratic functionals. Kerswell [41] demonstrated the duality of these approaches and showed that the optimal bounds obtained by these approaches exactly coincide.

Here, a new general method of obtaining rigorous bounds on various turbulent flow quantities is proposed, allowing the use of non-quadratic (but generally polynomial) functionals, and a preliminary estimate of its ability to improve the bounds obtained with quadratic functionals is made. The method is based on certain ideas that were introduced when the SOS polynomial analysis was applied to the stability problem; see, for example, [32] in the context of the analysis of global stability.

#### (i) The idea illustrated by a finite-dimensional approximation case

Many quantities of practical interest, for example the drag experienced by the body or the energy dissipation rate, are functionals of the velocity distribution. We denote such a quantity of interest with *Φ*, and we aim to find a bound for its time-averaged value.

Introducing a truncated Galerkin expansion
3.1
over an orthonormal solenoidal basis **u**_{i}(**x**) reduces the Navier–Stokes equations (2.1) to a system of ordinary differential equations
3.2
where
and the angular brackets denote a scalar product
Note that in this section the Galerkin expansion is used for the full velocity **u** and not for the velocity perturbation **u**′. The system (3.2) has a polynomial right-hand side, so it is natural to analyse its behaviour using the SOS techniques.

The total derivative of a function *V* (**a**) is defined as
Suppose that there exists some differentiable function *V* (**a**) and a constant *C* such that
3.3
We assume that the trajectories of the system remain bounded in some set, which is usually true for fluid flows, and we also require that the function *V* is bounded on the same set. Then, the average of d*V*/d*t* over infinitely long time is zero and, hence, the average of *Φ* over infinitely long time is bounded from below by *C*,
Consider the case when *Φ* is a polynomial function of *a*_{i}. For example, within approximation (3.1), the energy dissipation rate is
Then, one can search for a polynomial function *V* (**a**). In this case, *D*(**a**) is also a polynomial
3.4
As previously noted, verifying positive-definiteness of a polynomial is computationally very intensive. The computational complexity can be reduced significantly if the condition *D*(**a**)≥0 is replaced with a somewhat stronger condition that *D*∈*Σ*, where *Σ* is the set of all polynomials that can be represented as a sum of squares of other polynomials. Then, the problem of finding a good bound for can be formulated as an optimization problem over the coefficients of *V* (**a**) subject to the constraint *D*∈*Σ*,
3.5
For relatively modest polynomial degree and the number of variables, this problem can be solved using freely available Matlab toolboxes SOSTOOLS [20] or YALMIP [21,22].

#### (ii) Example

Consider now the nine-state example [34] modelling a Couette flow between two plates. We make the same assumptions as in [34,35]. For this system, the forcing term **F** is chosen in such a way that the steady state, that is, laminar, flow is the same for all values of the Reynolds number *Re*, but the forcing term depends on *Re*. Moreover, the basis is chosen in such a way that the steady flow is described by **u**_{1}. As a result, the steady flow is given by **a**=**h**=(1,0,…,0), and the force can be written accordingly as

Our calculations were carried out for the case of *L*_{x}=4*π* and *L*_{z}=2*π* of [35]. For this system, we solved the optimization problem (3.5) over a range of Reynolds numbers using either a second-degree function *V* in the form *V* =*p*_{2}(**a**) or a polynomial function in the form *V* =*p*_{2n}−*E*^{n+1}, where *E*=|**a**|^{2} and *p*_{n}(**a**) is a general *n*th-degree polynomial. The motivation for adding the term *E*^{n+1} is the same as in the stability analysis [32] and will not be discussed here. The bound *C** found by solving problem (3.5), where the optimization is performed over the coefficients of *p*_{2n}, is shown in figure 1 for *n* up to 3.

Also shown in figure 1 is the value of the steady-state dissipation rate *h*⋅*f*_{R}, which is also the upper bound for the dissipation rate in any flow, and the numerical estimate of computed by simulating the system (3.2) from random initial conditions and taking the time average of *Φ* over a suitably long interval. The numerical estimates use the average of 20 simulations at each Reynolds number.

While the background-flow method introduced in [40] appears to be rather different from the analysis here, in fact it is equivalent, up to the reduction to a finite-dimensional case, to selecting in our approach *V* =*p*_{2}. Because every positive-definite polynomial of a second degree can be represented as an SOS, our result for a second-degree polynomial *V* produces the best bound obtainable by the background-flow method applied to this finite-dimensional system. Comparing the results for the second-degree and higher-degree polynomials, one can judge the potential improvement in the bound owing to the change in degree. For *Re* below 40, the bound obtained with *V* =*p*_{4}−*E*^{3} practically coincides with the laminar values, whereas the bound obtained with the second-degree polynomial deviates by about 60% at *Re*≈40. However, for higher *Re*, the bounds become closer to each other. Increasing *n* from 2 to 3 gives very little improvement. This behaviour can be readily understood in view of the appearance of a periodic orbit at *Re*=80.54 (see, in particular, [35], fig. 7). This orbit is unstable and, therefore, it does not affect the simulated time-average of the dissipation rate. However, it does affect the bound. This is a generic problem for all attempts to approximate the actual turbulent dissipation rate with rigorous bounds: by their very nature, the bounds apply to both stable and unstable solutions. We suggest a way of dealing with this difficulty in §3*b*. The example given here shows that the approach based on the SOS might potentially produce much better bounds than the other known methods.

#### (iii) Extension to the full Navier–Stokes equation case

We start by projecting the Navier–Stokes equations (2.1) onto a solenoidal subspace, to eliminate the pressure. This gives
3.6
with subscript *S* denoting projection onto a solenoidal subspace. In practice, this projection would normally be achieved by Galerkin approximation.

Let *V* [**u**] be a functional of **u**. Then

Let *Φ*[**u**] also be a functional of **u**, for the time-average of which we would like to compute a bound. Denote

Because the time-averaged value of (d*V*/d*t*)[**u**] is zero, if *D*[**u**]≻0, then , and if then , where the bar denotes time-averaging. Taking *V* as a quadratic functional is very similar to the background-flow method of Doering & Constantin [40]. In general, one can obtain the best possible lower bound as a solution to the optimization problem *C* subject to *D*≻0. For the case of quadratic *V* [**u**], the solution of this optimization problem also generates the background flow giving the best possible bound. The change in the formulation needed to obtain the upper bound is obvious.

To solve the optimization problem using SOS, one can, in principle at least, attempt a direct approach similar to the approach described in §2*a*. We consider here, however, the approach based on the introduction of an uncertain system, as in §2*b*.

We introduce Galerkin approximation using the eigenfunctions of the energy stability problem as the basis and represent the flow velocity as a sum of *k* Galerkin modes and a residual field **u**_{s}(**x**), as given by (2.3). A minor difference is that in this case the expansion is done for the full flow velocity rather than for the perturbation velocity. Then, we proceed to obtain the uncertain system (2.6).

We then seek the functional *V* of the form *V* =*V* (*a*,*q*^{2}). Fairly straightforward, if somewhat lengthy, transformations give (compare with (2.7)):

Note that the correct choice of the expression for *Φ*[**u**] can be crucial. From energy conservation, it follows that the time average of the energy dissipation rate equals the time average of the work per unit time of the body force 〈**u**,**F**〉. Therefore, for the purpose of calculating averages over infinite time, we can replace the instantaneous energy dissipation rate with 〈**u**,**F**〉. It can then be represented as
Then,

Importantly, 〈**u**_{s},**F**〉 can be bounded via *q*^{2}, which would not be true for the instantaneous energy dissipation rate. Then, bounding also *Γ* and ** Θ** as in (2.6), we can formulate, similar to §2, an SOS problem in finite dimensions, looking for a polynomial

*V*(

*a*,

*q*

^{2}).

So far, this approach has not been applied to a particular flow. We return now to the finite-dimensional formulation in order to discuss more new ideas, not fully developed yet for the infinite-dimensional systems.

### (b) Eliminating the influence of unstable solutions on the bounds

The bounds obtained so far suffer from a serious limitation. Suppose that there is an unstable steady solution of (3.2): **a**(*t*)=**a**_{0}, **f**(**a**_{0})=0, where, for convenience, we rewrite (3.2) as
3.7
and rewrite (3.3) as
3.8
which implies that . An unstable solution does not affect the actual value of , because a small random noise, always present in reality, will result in the trajectory of the system leaving the vicinity of **a**_{0}. However, the existence of such an unstable solution does affect the lower bound given by (3.8), because, with **f**(**a**_{0})=0, no choice of *V* (**a**) can affect the value of *Φ*(**a**)−**f**(**a**)⋅∇_{a}*V* (**a**)−*C* at **a**=**a**_{0}. This problem is quite significant for obtaining bounds of parameters of turbulent flows, where the existence of an unstable laminar flow and/or unstable travelling wave can strongly affect the bound. Adding noise to (3.7) would resolve the problem, but it is not immediately obvious how one should modify (3.8) in this case. It turns out that progress can be made using a dual formulation.

#### (i) Dual formulation for the problem of bounds for time-averages

Assuming ergodicity, the time-average of *Φ* can be expressed as an ensemble-average over stationary probability density *ρ*(**a**),
where *ρ*≥0 satisfies the stationary Liouville equation
and the normalization condition
(The solution to this equation might turn out to be a generalized function, for example a delta function, but we ignore this now: with due caution, our arguments also apply in this case.) The lower bound on can be obtained as
3.9
The Lagrangian for this problem is
with λ(**a**)≥0. Assuming that *ρ*(**a**) decays fast enough as and using the Gauss theorem or integration by parts gives
The dual problem consists of maximizing
over *ν*, λ(**a**)≥0, and *μ*(**a**). Hence, for any *μ*(**a**) satisfying the regularity and behaviour-at-infinity constraints that ensure the applicability of the Gauss theorem, the expression
3.10
gives a lower bound on . Comparing this with (3.8) shows that our original formulation for the bound on the time-averaged quantity is dual to (3.9).

#### (i) Dual formulation for a system with random noise

Consider now the system
3.11
where the components *ξ*_{i}(*t*) are statistically independent delta-correlated Gaussian-distributed zero-mean random functions, with . Then, the steady probability density *ρ*(**a**) satisfies the steady Fokker–Planck equation
3.12
The Fokker–Planck equation replaces the Liouville equation as a constraint in (3.9),
3.13
Repeating the same arguments that were used in deriving (3.10) gives a lower bound for the system with noise as
3.14

Note that looking for a polynomial *μ*(**a**) using the SOS approach for the case with noise is formally not more complicated than doing this in the case of a zero noise, because adding the extra term does not increase the degree of the polynomial. We tested this approach with a rather trivial example. It is yet to be applied to a system of practical interest.

### (c) Flow control

As we have already mentioned in the Introduction, the SOS approach found extensive applications in control theory. Here, we concentrate specifically on a new idea for feedback control of turbulent flows. The final goal of such control is a reduction of turbulent drag or turbulent energy dissipation. For practically interesting regimes, however, fully suppressing turbulence, that is, making the laminar flow stable, is not realistic. This creates a difficulty for the mathematical approach to developing control laws for turbulent drag reduction, because, short of full direct numerical simulation of the turbulent flow, there is no way of calculating turbulent drag or dissipation. The existing works rely on semi-empirical, conceptual, correlations between turbulent energy dissipation and other features of the flow. This situation is best exemplified by the following citation from the well-known work of Lee *et al.* [42]: ‘We found that the choice of the cost functional to be minimized is critical in the performance of the control. Since the streamwise vortices have been known to be responsible for large drag in turbulent boundary layers, we tried to choose the cost functional that is directly related to them. This is based on our conjecture that a suitable manipulation of the streamwise vortices would lead to drag reduction’.

We propose to use the energy dissipation bound as the cost functional. This removes the empiricism and makes the entire scheme of determining the control law formal.

Just as an illustration, let us assume that the control system can be described by adding a polynomial function **G**(**a**), describing therefore the control law, to the right-hand side of (3.11),
If the lower bound in question is the lower bound for energy dissipation in turbulent flow driven by a given body force, then one can try to find a control aimed at increasing this bound in a hope that, when applied to the flow, this control also results in an increase in the energy dissipation rate itself (this is a beneficial change because, with the force fixed, this corresponds to a greater flow rate). With added control, **f**(**a**) should be replaced with **f**(**a**)+**G**(**a**) in (3.11). The bound is then
3.15
One can then try to design the best controller by optimizing over both *μ*(**a**) and **G**(**a**),

A known difficulty in applying the SOS approach in control is that, as one can see, *μ*(**a**) and **G**(**a**) taken together enter the cost functional nonlinearly. This issue, however, is well beyond the scope of this paper. Application of SOS in control is extensively described in the literature (e.g. [43–45]).

## 4. Conclusion

For finite-dimensional systems with a polynomial right-hand side, the SOS polynomial optimization approach allows one to construct Lyapunov functions in a formalized procedure carried out by a computer. Because the Navier–Stokes equations for incompressible flows have a quadratic nonlinearity, it is natural to attempt to apply SOS in fluid dynamics. The first part of this paper gives a review of the application of SOS to the problem of global stability of fluid flows. Direct application of this technique, described in §2*a*, requires derivation of an additional set of integral identities relating expressions involving functions and their partial derivatives. Doing this is not automated and remains dependent on the ingenuity of the researcher. An alternative approach is based on reducing the Navier–Stokes equations to an uncertain system, as described in §2*b*. Using this method, global stability of a version of the rotating Couette flow has already been studied, with the results providing a better stability range than the range given by the classic energy stability method. In this approach, an open question is how the quality of the available uncertainty bounds varies with the dimension of the uncertain system.

The second part of the paper contains new results and puts forward new ideas. In particular, a new method of obtaining bounds for time-averaged flow parameters, for example the energy dissipation rate, is proposed. The potential of the method was tested on the example of a model problem. It is shown that the proposed method of deriving a bound for a time-averaged quantity is a dual to the problem of obtaining a bound for the ensemble-averaged quantity. Adding random noise to the system and exploiting the duality, a method of obtaining bounds for time-averaged quantities, which is insensitive to unstable steady states and periodic orbits, is proposed. The method results in the SDP problem of the same level of difficulty as in the case of no noise. It is proposed to use the bound on the energy dissipation rate as the cost functional in the design of flow control aimed at reducing turbulent drag.

## Funding statement

Funding from EPSRC under grant no. EP/J011126/1 is gratefully acknowledged. D.H. is partially supported by the State Key Laboratory of Industrial Control Technology grant (ICT1304), Zhejiang University, PR China.

## Appendix A. Sum of squares methods and SOSTOOLS

Here, we provide more details on SOS methods and SOSTOOLS.

A polynomial *P*(*a*_{1},…,*a*_{n})≡*P*(**a**) is an SOS, if there exist polynomials *f*_{1}(**a**),…,*f*_{m}(**a**) such that
A 1
If there exists such a decomposition, then definitely *P*(**a**)≥0 for all . The converse is not always true, i.e. there are non-negative polynomials for which no SOS decomposition exists, but, in practice, replacing non-negativity with the existence of an SOS decomposition is adequate. However, testing whether *P*(**a**) is SOS is much more computationally tractable than testing non-negativity [1]: the former requires the solution of an SDP, whereas the latter is an NP-hard problem.

The way to test whether a decomposition of the form (A 1) exists, we write *P*(**a**) as explained in the main text,
where **m** is a vector of monomials in **a** of degree up to half the degree of *P*. For example, consider the polynomial
and therefore, by comparison, *q*_{11}=5,*q*_{22}=2,2*q*_{12}+*q*_{33}=−1,*q*_{13}=−1 and *q*_{23}=−1, subject to *Q*≥0: we have a semidefinite constraint on a matrix, subject to affine relationships in its entries, which is an SDP.

Often, we want to establish ‘conditional satisfiability’ conditions of the form
To address such questions, we can use existing results from real algebraic geometry (*Positivstellensatz*) and search for SOS multipliers such that
This is true, as the summation is non-positive when *p*_{j}(**a**)≤0; hence, in that case, *p*_{0}(**a**)≥0. Such constraints can also be represented as above, as they encompass positive semidefinite constraints on unknown matrices with affine relationships on their entries.

Software exists to formulate SOS constraints as SDPs (SOSTOOLS [20] and YALMIP [21]), which can then be solved using SDP solvers, e.g. LMILAB, Mosek, PENBMI and PENSDP, CSDP, KYPD, SDPA, 8, SeDuMi, SDPLR and SDPNAL. SOS programs of the most general form that can be handled by SOSTOOLS take the form:

*Optimization*: minimize the linear objective function
where **c** is a vector formed from the (unknown) coefficients of polynomials *q*_{i}(**a**), for and SOS *q*_{i}(**a**), for which need to be found such that
In the above formulation, **w** is the vector of weighting coefficients for the linear objective function.

To define and solve an SOS program (SOSP) using SOSTOOLS, the steps include

Initialize the SOSP.

Declare the SOSP variables.

Define the SOSP constraints.

Set objective function (for optimization problems).

Call solver.

Obtain solutions.

Please check the user manual for a range of examples and customized functions.

## Footnotes

One contribution of 15 to a Theme Issue ‘Stability, separation and close body interactions’.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.