Sum-of-squares approach to feedback control of laminar wake flows

Abstract

In this paper a novel nonlinear feedback control design methodology for incompressible fluid flows aiming at the optimisation of long-time averages of flow quantities is presented. It applies to reduced-order finite-dimensional models of fluid flows, expressed as a set of first-order nonlinear ordinary differential equations with the right-hand side being a polynomial function in the state variables and in the controls. The key idea, first discussed in Chernyshenko et al. (Phil. Trans. R. Soc. Lond. A, vol. 372, 2014, 20130350), is that the difficulties of treating and optimising long-time averages of a cost are relaxed by using the upper/lower bounds of such averages as the objective function. In this setting, control design reduces to finding a feedback controller that optimises the bound, subject to a polynomial inequality constraint involving the cost function, the nonlinear system, the controller itself and a tunable polynomial function. A numerically tractable and efficient approach to the solution of such optimisation problems, based on sum-of-squares techniques and semidefinite programming, is proposed. To showcase the methodology, the mitigation of the fluctuation kinetic energy in the unsteady wake behind a circular cylinder in the laminar regime at $Re=100$, via controlled angular motions of the surface, is numerically investigated. A compact reduced-order model that resolves the long-term behaviour of the fluid flow and the effects of actuation, is first derived using proper orthogonal decomposition and Galerkin projection. In a full-information setting, feedback controllers are then designed to reduce the long-time average of the resolved kinetic energy associated with the limit cycle. These controllers are then implemented in direct numerical simulations of the actuated flow. Control performance, total energy efficiency and the physical control mechanisms identified are analysed in detail. Key elements of the methodology, implications and future work are finally discussed.

1 Introduction

In the past decades, the coordinated efforts of laboratory experiments using high-resolution flow diagnostics and large-scale direct numerical simulations (DNS) have considerably progressed our understanding of key physical processes and mechanisms in turbulent flows. Despite these new discoveries, progress in the ability to effectively control their spatiotemporal evolution in complex geometries has remained more elusive, owing to the nonlinear, multiscale nature of turbulent motion. Interest in control is driven by the huge economic, societal and environmental benefits that advances in the field would provide. Hence, the development of active control strategies is commonly regarded as one of the key enablers for future advances in efficient transportation, energy generation and distribution, and in many other technologically relevant industrial sectors.

Controlling and mitigating large-scale velocity fluctuations in the flow around bluff bodies, the problem that we discuss in this paper, is one such instance. When the Reynolds number exceeds a critical value, the periodic generation and shedding of organised vortical structures from the body produces intense fluctuations in the aerodynamic forces, resulting in structural fatigue (Sarpkaya 2004), acoustic noise production (Blevins 1984; Thomas, Kozlov & Corke 2008) and other undesirable effects, such as vortex-induced vibrations (Williamson & Govardhan 2004). The technological relevance of these flows has thus spawned significant interest in devising control strategies to tame their evolution. A variety of actuation/sensing strategies and control design methods have been proposed, as recently reviewed by Choi, Jeon & Kim (2008).

Because of the simplicity of the geometry, the two-dimensional flow past a circular cylinder has become the paradigmatic flow model to investigate vortex dynamics around bluff bodies. The laminar, steady solution is characterised by two recirculation eddies, whose length grows linearly with $Re$ (Fornberg 1985) and becomes unstable in a Hopf bifurcation at $Re\approx 47$ (Provansal, Mathis & Boyer 1987; Noack & Eckelmann 1994) due to a symmetry-breaking unstable global mode (Tang & Aubry 1997). The ensuing nonlinear regime saturates in a sustained periodic motion, vortex shedding, a stable periodic orbit attracting trajectories in an appropriately defined phase space of the system (Rempfer 2000; Noack et al. 2003) before the occurrence of other bifurcations at higher $Re$ (Barkley & Henderson 1996).

Control of this specific regime became a useful benchmark problem to develop and test novel feedback control design strategies (Lehmann et al. 2005). One of the perspectives adopted in several investigations on control has been the stabilisation by feedback of the unstable, steady, laminar wake flow. At low supercritical Reynolds numbers, only one unstable global mode, the Kármán mode, exists. Hence, proportional control strategies, where the signal from a single sensor located at some point in the wake is multiplied by a constant gain and fed back to the actuator, have been considered extensively (e.g. Berger 1967; Monkewitz, Berger & Schumm 1991; Roussopoulos 1993; Park, Ladd & Hendricks 1994). In the light of DNS and reduced-order modelling techniques for linear systems, Illingworth, Naito & Fukagata (2014) review succinctly some of these efforts and discuss the ‘gain window effect’ observed in previous numerical and experimental works, i.e. when suppression of the wake instability is achieved only if the gain is within a certain interval. They show that such an effect does not result from the destabilisation by control of other unstable modes, but rather that it is driven by the properties of the closed-loop system, in particular by the time delays associated with the feedback arrangement. The authors also showed that the window shrinks as $Re$ increases and it does not exist any more at $Re=80$ , reflecting the objective difficulty or impossibility of obtaining stabilising controllers as the Reynolds number increases. They concluded by suggesting that better control strategies, with more complicated dynamics than proportional control, might be required to improve performance.

Camarri & Iollo (2010) proposed a linear feedback design method, for the flow past a square cylinder, based on linearised dynamics and global linear stability analysis of the equilibrium solution, inspired by previous works on passive control design methods (see Giannetti & Luchini 2007; Marquet, Sipp & Jacquin 2008, and references therein). Camarri & Iollo (2010) proceed by examining the sensitivity of the linear stability problem with respect to the controller parameters, in order to displace the eigenvalues of the unstable and least stable modes to the left half of the complex plane via control design. They pointed out that the performance of this controller far from the design state, i.e. the control of the nonlinear saturated regime, needs to be explored a posteriori. They show that their feedback strategy can stabilise the fully nonlinear regime up to twice the critical Reynolds number of the natural flow. At higher Reynolds number, in highly nonlinear regimes, performance worsens. Interestingly, the authors point out that the basin of attraction of the stabilised wake structure shrinks consistently as the Reynolds number increases.

Carini, Pralits & Luchini (2015) investigated feedback control in the framework of linear optimal control theory, and designed and tested a full-dimensional minimum-control-energy compensator, free from spillover effects induced by the excitation by actuation of stable dynamics, often observed when control is designed on a reduced-order model (Barbagallo, Sipp & Schmid 2009). Using the feedback from a single sensor measuring the cross-stream velocity to control the rotation rate of the cylinder around its axis, Carini et al. (2015) showed that complete stabilisation of the unstable mode was possible only up to $Re\approx 59$ , if the sensor was located in a narrow region between 2 and 2.5 diameters downstream of the cylinder axis. The critical $Re$ was increased to 72 when a full-information controller was employed. They commented on this difference by suggesting that better performance on the nonlinear saturated flow could be obtained by adopting a nonlinear observer and ultimately a nonlinear control strategy.

These and other investigations have demonstrated that, as the Reynolds number increases, the flow dynamics becomes so strongly nonlinear to render linearisation of the equations around the unstable equilibrium and linear design methods scarcely effective. In the terminology of Brunton & Noack (2015), such systems are ‘irreducible’, in the sense that key nonlinear processes, such as vortex pairing/merging, inter-modal energy transfers and advection of coherent structures, crucial to describe the developed state of natural instabilities that arise progressively as the Reynolds number increases, cannot be described by a linearised theory. Furthermore, the gradual loss of linear stabilisability as the Reynolds number increases (Lauga & Bewley 2003), coupled with sensing/actuation constraints of practical technological nature, suggest that the developed nonlinear state of the flow needs to be addressed directly in the design stage.

Strategies where the structure of the feedback controller is heuristically fixed a priori and appropriate gains are obtained from optimisation or parameter exploration over nonlinear controlled regimes have been proposed (e.g. Fujisawa, Kawaji & Ikemoto 2001; Siegel, Cohen & McLaughlin 2006; Weller, Camarri & Iollo 2009; Lu et al. 2011; Mehmood et al. 2014). Weller et al. (2009) introduced a feedback structure consisting of a linear proportional controller relating several cross-flow velocity measurements in the near wake to the signal driving the actuators, two blowing/suction slots on the top and bottom walls of the square cylinder arrangement driven in opposite phase. Optimisation of the gains, to reduce the short-time-averaged $L_{2}$ -norm of the difference between the instantaneous flow field and the unstable steady solution, was then performed in a trust-region, reduced-order, adaptive setting. The resulting feedback arrangement was able to stabilise the flow starting from the saturated nonlinear regime at a Reynolds number almost twice the critical value. However, because the optimisation involved a cost function defined over a finite short horizon, the best controller resulted in excellent performance in this interval but performance subsequently degraded, especially at larger Reynolds numbers. Weller et al. (2009) concluded by pointing out that the asymptotic stability of the closed-loop system cannot be ensured by their method, as the long-term behaviour of the system is not considered in the design.

These investigations strongly relied on the ingenuity of the researchers, on heuristic choices of sensor/actuation position and type, and on solid understanding of the flow physics. Such heuristic strategies might show significant limitations when applied to flows with more complex nonlinear dynamics. Recent model-free approaches, such as genetic programming control (see e.g. Debien et al. 2016; Parezanović et al. 2016, and references therein), use evolutionary strategies to automatically discover such heuristics in experimental control studies, using a black-box optimisation approach. These approaches can lead to the emergence of unexpected control solutions, as they effectively explore large search spaces, and can uncover novel control mechanisms (Gautier et al. 2015).

On the other side of the spectrum, optimal control theory (Abergel & Temam 1990) is probably one of the most versatile model-based control design methods for nonlinear systems. Optimal control, in the predictive setting, consists in finding and applying in a feed-forward fashion the control input that optimises a suitable cost function defined as a definite integral spanning a predetermined finite horizon. Although such a strategy is extremely computationally expensive, it is considered to represent the upper limit on the achievable control performance (Bewley, Moin & Temam 2001). Optimal control of a circular cylinder wake via rotary actuation has been implemented by Protas & Styczek (2002) to minimise a cost function involving the sum of the power associated with control and that associated with the drag, using optimisation horizons up to roughly one vortex shedding period. More recently, Flinois & Colonius (2015) implemented the same algorithm but significantly extended the optimisation horizon, up to 100 convective time units, i.e. at least 10 times larger than previous efforts. The important observation is that long-time horizons, representative of the long-term behaviour of the controlled system, were necessary to suppress vortex shedding at Reynolds numbers between 75 and 200, and achieved far better performance, with smoother control inputs, than the short-time horizon approach of Protas & Styczek (2002), enabling the identification of physical control mechanisms.

1.1 Objectives and structure of the paper

The main purpose of this paper is to present a novel paradigm for model-based feedback control of fluid flows, in an effort to address some of the outstanding issues discussed in the introduction. Firstly, the proposed control paradigm applies directly to nonlinear Galerkin-type models of incompressible fluid flows. No linearisation around an operating point is performed and the only dynamical approximation is the truncation of the Galerkin velocity expansion. Hence, important nonlinear processes that can be described by such models can be controlled, if not exploited. Secondly, the long-term behaviour of the system is central in the design, as the optimisation targets long-time averages, defined over infinite horizons. The key step to overcome the objective difficulty of treating and optimising such averages is to replace it by estimation/optimisation of bounds, as first proposed by Chernyshenko et al. (2014).

The theoretical and algorithmic backbone of this approach is a recent breakthrough in control theory and optimisation, i.e. the discovery that the sum-of-squares (SOS) decomposition of a polynomial can be found, if it exists, via the solution of a semidefinite program (SDP) (Parrilo 2003). These advances have recently emerged as a promising basis to solve many computationally hard analysis/design problems for systems whose dynamics are described by polynomial functions, such as the estimation of the attraction region of equilibria (Valmorbida, Tarbouriech & Garcia 2013) as well as the simultaneous optimisation of a stabilising controller and a high-degree Lyapunov function certifying the stability of the controlled system (Zhao & Wang 2010; Nguang, Krug & Saat 2011). These new paradigms of design and analysis provide us with numerically tractable methods to take a new perspective of many fundamental problems in fluid dynamics as reviewed in Chernyshenko et al. (2014), such as nonlinear control design, the objective of this paper, or nonlinear stability analysis, as in Huang et al. (2015).

The paper is organised as follows. In § 2, a concise presentation of SOS techniques is reported. Numerous references to this technology are presented for the more interested reader. Section 3 describes the control design methodology, using a relatively general notation. More specifically, it discusses the technique used to estimate bounds on long-time averages and its application to control design via bounds optimisation. In § 4 these ideas are applied to the benchmark control problem of regulation of vortex shedding past a circular cylinder at a Reynolds number equal to 100, via rotary oscillations of the surface. This problem was extensively discussed in this introduction to put our results in a more focused context and was chosen as a pretext to describe a methodology that applies independently of the specific case, i.e. from the details of the flow, the actuation/sensing arrangement and the modelling approach. The numerical set-up is discussed first. The model order reduction strategy, based on proper orthogonal decomposition (POD) and Galerkin projection, is then introduced. State-feedback controllers are further designed and performance is assessed by implementation in DNS in a full-information setting. Conclusions and future work to be addressed are summarised in § 6.

2 The sum-of-squares decomposition

We provide in this section a succinct overview of SOS techniques, in order to convey the general underlying ideas. In this section, we favour clarity over mathematical rigour, with the hope of bridging the gap between the mathematical aspects and fluid mechanics. We refer the interested reader to our previous work (Chernyshenko et al. 2014), where a broader overview of the SOS technique and its applications in fluid mechanics is given.

Despite the complexity of the underlying mathematical framework, the idea of the SOS decomposition of a polynomial is rather simple. As an example, one might be interested in checking the non-negativity of a given multivariate polynomial function $P(a_{1},\ldots ,a_{N})=P(\boldsymbol{a})$ , of even degree $2d_{P}$ , that is, if $P(\boldsymbol{a})\geqslant 0$ for all $\boldsymbol{a}\in \mathbb{R}^{N}$ . Checking non-negativity for a general multivariate polynomial is NP-hard, hence intractable from a computational perspective (Papachristodoulou & Prajna 2002). However, a sufficient condition for $P(\boldsymbol{a})$ to be non-negative is that it can be decomposed into the sum of the squares of $M$ polynomial functions $p_{1}(\boldsymbol{a}),\ldots ,p_{M}(\boldsymbol{a})$ , of lower degrees as

(2.1) $$\begin{eqnarray}P(\boldsymbol{a})=\mathop{\sum }_{i=1}^{M}p_{i}^{2}(\boldsymbol{a}).\end{eqnarray}$$

Finding such a decomposition is equivalent to finding a positive semidefinite matrix $\unicode[STIX]{x1D64D}$ , which can be assumed symmetric without loss of generality, and a suitable vector $\boldsymbol{v}(\boldsymbol{a})$ containing monomials in the variables $a_{i}$ up to degree $d_{P}$ such that

(2.2) $$\begin{eqnarray}P(\boldsymbol{a})=\boldsymbol{v}(\boldsymbol{a})^{\text{T}}\unicode[STIX]{x1D64D}\boldsymbol{v}(\boldsymbol{a}).\end{eqnarray}$$

If one can find a positive semidefinite $\unicode[STIX]{x1D64D}$ , then a linear transformation of coordinates can reduce it to a diagonal form, with non-negative entries on the main diagonal, reducing $P$ to a linear combination of squares of polynomials, clearly being equivalent to non-negativity. However, the converse is not necessarily true, that is, not all non-negative polynomials admit an SOS decomposition, a famous counter-example being the Motzkin polynomial.

In a design problem, it might be of interest to construct a non-negative polynomial function subject to a set of constraints, rather than checking non-negativity of an existing one. This problem, which we will deal with in what follows, can be treated essentially using the same approach. It is worth noting that, in practice, the decomposition (2.2) is only approximate, and the error

(2.3) $$\begin{eqnarray}e(\boldsymbol{a})=P(\boldsymbol{a})-\boldsymbol{v}(\boldsymbol{a})^{\text{T}}\unicode[STIX]{x1D64D}\boldsymbol{v}(\boldsymbol{a})\end{eqnarray}$$

is non-zero. However, by Theorem 4 in Löfberg (2009), the polynomial $P(\boldsymbol{a})$ is still certifiably non-negative if

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{min}-\dim (\unicode[STIX]{x1D64D})\times |r|\geqslant 0,\end{eqnarray}$$

where $\unicode[STIX]{x1D706}_{min}$ is the smallest eigenvalue of $\unicode[STIX]{x1D64D}$ , $\dim (\unicode[STIX]{x1D64D})$ denotes the dimension of the matrix $\unicode[STIX]{x1D64D}$ , and $r$ is the coefficient of $e(\boldsymbol{a})$ that has the largest magnitude.

From a computational perspective, finding the SOS decomposition of a polynomial amounts to finding a positive semidefinite matrix $\unicode[STIX]{x1D64D}$ , subject to a set of linear equality constraints, arising from the equality in (2.2). This problem can be reformulated as an SDP (Parrilo 2003), a convex and tractable problem to solve. Several freely available software tools that can formulate and solve efficiently this kind of problem exist, such as the Matlab toolboxes SOSTOOLS (Papachristodoulou et al. 2013) and YALMIP (Löfberg 2004).

3 The control design method

3.1 Problem statement

We consider finite-dimensional dynamical systems given as a set of nonlinear coupled ordinary differential equations, as

(3.1) $$\begin{eqnarray}\frac{\text{d}\boldsymbol{a}}{\text{d}t}=\,\boldsymbol{f}(\boldsymbol{a},u),\end{eqnarray}$$

where $\boldsymbol{a}\in \mathbb{R}^{N}$ is the state-variables vector, $u\in \mathbb{R}$ is the control, and $\boldsymbol{f}:\mathbb{R}^{N}\times \mathbb{R}\rightarrow \mathbb{R}^{N}$ is assumed to be a polynomial function in the state variables and in the control. For the sake of reducing clutter, we discuss a single input case, but the multiple input case can be treated with minor revisions in the derivation. For incompressible fluid flows, this is the formulation that results naturally from Galerkin projection of the governing equations onto a finite-dimensional orthonormal set of basis functions (Fletcher 1984). It is well known that, for such systems, the vector field $\boldsymbol{f}$ can have constant, linear and quadratic terms, and the latter conserves energy for a large class of boundary conditions of the original partial differential equations. The way the control appears in $\boldsymbol{f}$ depends on the type of actuation: for volume forces, the right-hand side is affine with the input $u$ ; whereas for actuation via the boundary, a ‘lifting’ procedure results in the dynamics of the system being dependent on $\text{d}u/\text{d}t$ and $u^{2}$ .

Suppose that for system (3.1) it is of interest to control the value of a turbulent quantity $\unicode[STIX]{x1D6F7}(t)$ , the cost. This could express, for instance, the instantaneous turbulent kinetic energy, or the energy dissipation rate. Suppose further that the cost can be expressed as a positive-definite polynomial function of the state variables and of the control, i.e. $\unicode[STIX]{x1D6F7}(t)=\unicode[STIX]{x1D6F7}(\boldsymbol{a}(t),u(t))$ . In general, but more specifically for systems exhibiting turbulent behaviour, long-time statistics of $\unicode[STIX]{x1D6F7}(t)$ , for example, long-time averages,

(3.2) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}}=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\unicode[STIX]{x1D6F7}(\boldsymbol{a}(t),u(t))\,\text{d}t,\end{eqnarray}$$

are of primary interest, where $\boldsymbol{a}(t)$ is the solution of (3.1) with $u=u(t)$ and for some initial condition. Denoting first by $\overline{\unicode[STIX]{x1D6F7}}^{0}$ the long-time-averaged cost without control, the objective is to design a state-feedback controller

(3.3) $$\begin{eqnarray}u(t)=g(\boldsymbol{a}(t))\end{eqnarray}$$

that manipulates the long-term behaviour of (3.1) such as to reduce, or increase depending on the problem, the long-time-averaged cost to $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ . Here, we also restrict $g:\mathbb{R}^{N}\rightarrow \mathbb{R}$ to be an initially undetermined polynomial function of arbitrary degree $d_{g}$ in the state variables, in order to leverage SOS techniques. Such a restriction imposes a high degree of smoothness on the control, but highly nonlinear controllers can be designed, as $d_{g}$ can be regarded as a design parameter. Note that the controller (3.3) makes the closed-loop system (3.1) an autonomous system. We assume here that complete and exact information on the instantaneous state of the system is available; hence, we avoid the necessity of designing an observer. This step would be required in a practical application, but it is beyond the scope of this paper, which focuses on control design only.

Figure 1. Illustration of the general idea behind the proposed control methodology. Instead of designing a controller that reduces the time average from $\overline{\unicode[STIX]{x1D6F7}}^{0}$ to $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ , a controller that reduces the upper bound from $C^{0}$ to $C^{\ast }$ is sought. Under the action of such a controller, the time average is also expected to decrease, although this cannot be guaranteed in a general case.

Ideally, the controller could be designed by solving the optimisation problem

(3.4) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}}^{\ast }=\left\{\begin{array}{@{}l@{}}\displaystyle \min _{g}~\overline{\unicode[STIX]{x1D6F7}}\quad \\ \text{subject to}~\dot{\boldsymbol{a}}=\boldsymbol{f}(\boldsymbol{a},g(\boldsymbol{a})),\quad \end{array}\right.\end{eqnarray}$$

where $\dot{\boldsymbol{a}}=\text{d}\boldsymbol{a}/\text{d}t$ . The non-convexity of (3.4), but most importantly the fact that the minimisation of long-time averages are considered, makes (3.4) difficult to solve. The key step, previously suggested in Chernyshenko et al. (2014), is illustrated in figure 1. Instead of treating a long-time average directly, we replace the original problem with the analysis of an upper bound of the average, i.e. a value $C$ for which an algorithm exists proving $\overline{\unicode[STIX]{x1D6F7}}\leqslant C$ for system (3.1), where the equality holds when the bound is tight. Hence, instead of attempting to reduce the long-time average, we reformulate (3.4) into the problem of designing a controller minimising the upper bound, from $C^{0}$ , the bound on $\overline{\unicode[STIX]{x1D6F7}}^{0}$ , to $C^{\ast }$ , the bound on $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ . This reads as

(3.5) $$\begin{eqnarray}C^{\ast }=\left\{\begin{array}{@{}l@{}}\displaystyle \min _{g}~C\quad \\ \text{subject to}~\bar{\unicode[STIX]{x1D6F7}}\leqslant C,~\dot{\boldsymbol{a}}=\boldsymbol{f}(\boldsymbol{a},g(\boldsymbol{a})).\quad \end{array}\right.\end{eqnarray}$$

The hope is that, under the action of such a controller, the actual time average $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ will also decrease. This is not guaranteed to happen in a general case. As a trivial, yet representative, example, consider a system having multiple stable equilibria $\boldsymbol{a}_{i}$ , each with its own basin of attraction. In such a case, the long-term behaviour of trajectories, hence the time average, depends on the particular choice of the initial conditions. An upper bound on the time average of some cost $\unicode[STIX]{x1D6F7}(\boldsymbol{a})$ is simply $\max _{i}\unicode[STIX]{x1D6F7}(\boldsymbol{a}_{i})$ . The crucial point is that a controller designed to reduce the upper bound is guaranteed to decrease the actual time average only in a ‘worst-case scenario’, i.e. when the trajectory starts from the basin of attraction of the steady solution associated with the bound. Should this not be the case, it is perfectly possible that the actual long-time average will increase.

The occurrence of such a behaviour depends on the particular choice of the cost function $\unicode[STIX]{x1D6F7}(\boldsymbol{a})$ and on the topology of the system’s phase space, i.e. the attractors and/or repellors that populate it. Nevertheless, manipulating and analysing bounds is much easier than doing so with long-time averages directly. SOS techniques can be employed exactly for such a purpose, as we show in the next section. In the case when the algorithm used to calculate the upper bound does not guarantee that the bound is tight, the outcome of the optimisation depends also on the algorithm, which, therefore, should be specified. This is done in the following section.

3.2 Bounds estimation

The first step is to derive an upper bound $C^{0}$ on the average $\overline{\unicode[STIX]{x1D6F7}}^{0}$ , for uncontrolled dynamics. We define a polynomial function in the state variables, $V(\boldsymbol{a})$ , of degree $d_{V}$ , containing unknown decision variables as its coefficients. We assume that trajectories of the system (3.1) are bounded in some set, as one would expect in a dissipative system such as a fluid flow, both in the infinite-dimensional case, and for non-degenerate finite-dimensional representations thereof.

The function $V$ is also bounded in this set, as it is a polynomial. The total time derivative of $V$ along trajectories of the system,

(3.6) $$\begin{eqnarray}\frac{\text{d}V(\boldsymbol{a})}{\text{d}t}=\frac{\unicode[STIX]{x2202}V}{\unicode[STIX]{x2202}\boldsymbol{a}}\boldsymbol{\cdot }\frac{\text{d}\boldsymbol{a}}{\text{d}t}=\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a}),\end{eqnarray}$$

is then also bounded, where $\unicode[STIX]{x1D735}_{\boldsymbol{a}}V\triangleq \unicode[STIX]{x2202}V/\unicode[STIX]{x2202}\boldsymbol{a}$ is the gradient of $V$ with respect to the coordinates of the phase space. Now, suppose one can find some $V$ such that the following polynomial inequality

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a})+\unicode[STIX]{x1D6F7}(\boldsymbol{a})\leqslant C\end{eqnarray}$$

is satisfied for all $\boldsymbol{a}$ in $\mathbb{R}^{N}$ , for a given $C$ . Then, it is straightforward to show that $C$ is an upper bound for $\overline{\unicode[STIX]{x1D6F7}}^{0}$ , i.e. $\overline{\unicode[STIX]{x1D6F7}}^{0}\leqslant C$ . This is because, when taking the time average of (3.7) with $\boldsymbol{a}=\boldsymbol{a}(t)$ , the term

(3.8) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a})}=\overline{\,\frac{\text{d}V(\boldsymbol{a})}{\text{d}t}\,}=\lim _{T\rightarrow \infty }\frac{1}{T}\int _{0}^{T}\frac{\text{d}V(\boldsymbol{a})}{\text{d}t}\text{d}t=\lim _{T\rightarrow \infty }\frac{1}{T}[V(\boldsymbol{a}(T))-V(\boldsymbol{a}(0))]\end{eqnarray}$$

vanishes identically under the above assumption of boundedness. Hence, the upper bound $C^{0}$ can be obtained by minimising $C$ over all possible polynomials $V$ of a given degree under the polynomial constraint (3.7), i.e. by solving

(3.9) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}}^{0}\leqslant C^{0}=\left\{\begin{array}{@{}l@{}}\displaystyle \min _{V}~C\quad \\ \text{subject to}~-(\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a})+\unicode[STIX]{x1D6F7}(\boldsymbol{a})-C)\geqslant 0.\quad \end{array}\right.\end{eqnarray}$$

Because verifying positive-definiteness of a given polynomial, as well as constructing one as in the present case, is a notoriously difficult problem, we replace the constraint in (3.9) such as to have

(3.10) $$\begin{eqnarray}\overline{\unicode[STIX]{x1D6F7}}^{0}\leqslant C^{0}\leqslant C_{SOS}^{0}=\left\{\begin{array}{@{}l@{}}\displaystyle \min _{V}~C\quad \\ \text{subject to}~-(\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a})+\unicode[STIX]{x1D6F7}(\boldsymbol{a})-C)\in \unicode[STIX]{x1D6F4},\quad \end{array}\right.\end{eqnarray}$$

where $\unicode[STIX]{x1D6F4}$ is the set of all polynomials that have an SOS decomposition. From a numerical point of view, this optimisation problem is solved by trial and error by reducing $C$ until a $V$ satisfying the polynomial inequality cannot be found. For a given $C$ , the search for the function $V$ is numerically reformulated into an SDP using standard software tools (Löfberg 2004; Papachristodoulou et al. 2013). It is a convex problem, hence can be solved efficiently, and the solution, if it exists, is unique. In general, a hierarchy of bounds can be obtained by increasing the degree of the polynomial function $V$ . Note that the same procedure can be used to estimate a lower bound, when maximisation of the time average is of interest, by reversing the sign of the inequality in (3.7), and change (3.10) to a maximisation problem.

Strengthening the non-negativity constraint to an SOS constraint adds conservativeness in the optimisation, in the sense that the upper bound found from (3.10) can be, in principle, lower than the bound that could be found if one was able to solve (3.9) directly, so the tightness of the obtained bound may not be guaranteed. This is because not all positive-definite polynomials can be decomposed into the sum of squares of other polynomials, although this appears to be a special case (Tan 2006). However, the second problem is numerically tractable, whereas the former is not.

It is worth pointing out that finding a finite upper bound of a long-time average on a positive-definite cost, with the method defined above, does not automatically prove the boundedness of the trajectory of the system. Bounds on long-time averages are determined by the invariant sets that populate the phase space of the system, no matter what their stability is, because they define the permanent regime. The bound could be given by an unstable set, e.g. a repellor, and, in the absence of other information, it is not possible to assert boundedness of the trajectories. A modification of the polynomial inequality (3.9), based on the idea of adding stochastic noise to (3.1), to include only stable invariant sets, is discussed in Chernyshenko et al. (2014). However, with this modification, a finite upper bound only proves boundedness of trajectories for almost all initial conditions. In fact, there might still exist an unbounded unstable set, which might allow a trajectory with a particular initial condition, on this set, to escape to infinity. The inability to find an upper bound does not prove that trajectories are unbounded. This is because the SOS constraint in (3.10) is stronger than the non-negativity constraint in (3.9), resulting in $C^{0}\leqslant C_{SOS}^{0}$ . Hence, one could probably formulate a corner-case problem where a finite $C_{SOS}^{0}$ cannot be found, whereas $C^{0}$ exists and is finite.

A recent discussion on such issues is given by Schlegel & Noack (2015). These authors proposed a computational procedure that can be used to prove boundedness of the trajectories. It is based on the idea of finding a globally attracting ‘trapping region’, i.e. a closed set in the phase space such that all trajectories converge to this region and remain inside it once they have entered it. Their procedure is based on finding an appropriate shift in the phase space, such that the perturbation energy in this translated reference frame possesses the mathematical properties of a Lyapunov function for large deviations from the shifted origin. From a computational perspective, they employed a simulated annealing algorithm, or ad hoc searches along particular directions in the phase space, to identify the appropriate shift. However, this algorithm can only prove the existence of a trapping region – it cannot disprove it. In appendix A we report an alternative and rigorous SOS-based procedure that can be used to prove the existence of a monotonically trapping region.

3.3 Bounds optimisation

As for the bound estimation problem, we consider a tunable polynomial function $V(\boldsymbol{a})$ , and assume initially that trajectories remain bounded under closed-loop control. The optimisation problem equivalent to (3.10) is now

(3.11) $$\begin{eqnarray}C^{\ast }\leqslant C_{SOS}^{\ast }=\left\{\begin{array}{@{}l@{}}\displaystyle \min _{V,g}~C\quad \\ \text{subject to}~-(\unicode[STIX]{x1D735}_{\boldsymbol{a}}V(\boldsymbol{a})\boldsymbol{\cdot }\boldsymbol{f}(\boldsymbol{a},g(\boldsymbol{a}))+\unicode[STIX]{x1D6F7}(\boldsymbol{a})-C)\in \unicode[STIX]{x1D6F4},\quad \end{array}\right.\end{eqnarray}$$

where the minimisation of the upper bound is now performed over all possible polynomial functions $V(\boldsymbol{a})$ and state-feedback polynomial controllers $g(\boldsymbol{a})$ , of given degree $d_{V}$ and $d_{g}$ , respectively. The additional degrees of freedom associated with $g$ can allow, unless one is dealing with certain pathological cases, a further reduction of the upper bound, that is, $C_{SOS}^{\ast }<C_{SOS}^{0}$ . As previously anticipated, it is not guaranteed that the feedback controller obtained using this procedure will reduce the actual value of the time average of the cost in closed-loop simulation of the system, that is, the inequality $\overline{\unicode[STIX]{x1D6F7}}^{\ast }<\overline{\unicode[STIX]{x1D6F7}}^{0}$ cannot be guaranteed to hold.

The bounds $C^{0}$ and $C^{\ast }$ solely depend on the analytic definition of the vector field $\boldsymbol{f}$ , hence on the structure of the system’s phase space. Because the system’s invariant sets determine the long-term evolution, hence the bound, one can see this design scheme as finding the vector field induced by $g(\boldsymbol{a})$ that moves/reshapes the set associated with the bound such as to reduce favourably the long-time-averaged cost.

The bound optimisation problem is still non-convex, because one needs to optimise simultaneously the tunable function $V(\boldsymbol{a})$ and the controller $g(\boldsymbol{a})$ , and so the tuning variables in $V$ are multiplied by those in $g$ . This problem is not directly reducible to an SDP and convex optimisation techniques cannot be readily applied. This is a well-known problem in the SOS community, and is similar to that encountered when optimising simultaneously a globally stabilising feedback controller and a polynomial Lyapunov function certifying the global stability of an equilibrium (Zhao & Wang 2010; Nguang et al. 2011). Alternative iterative algorithms need to be used (see e.g. Henrion & Garulli 2005, for an overview).

In this paper we have developed a similar iterative algorithm, which is described in detail in appendix B and used to solve (3.11). The details are as follows.

1. (1) For a given upper bound $C$ and an initial controller $g(\boldsymbol{a})$ , whose derivative can be calculated in a certain way, e.g. equation (5.1) when $g(\boldsymbol{a})$ is linear, check the feasibility of the resultant SOS problem, namely, equation (2.4), by tuning $V(\boldsymbol{a})$ . Here, the feasibility of SOS optimisation implies that the controller $g(\boldsymbol{a})$ is effective in reducing the upper bound to $C$ .

2. (2) Fixing the optimised $V(\boldsymbol{a})$ and still keeping $\text{d}g/\text{d}t$ as in (1), further minimise the upper bound $C$ by solving the resulting SDPs in the decision variables of $g(\boldsymbol{a})$ . Note that the Schur complement technique will be adopted to resolve the nonlinearity of the SOS problem, as demonstrated in (B 3). In addition, only a reduction of $\unicode[STIX]{x1D6FF}C$ is considered at each step.

3. (3) Update $\text{d}g/\text{d}t$ using the optimised $g(\boldsymbol{a})$ in (2) and repeat the procedure (1)–(2) until $C$ cannot be decreased any more.

4. (4) Output the optimised $C$ and the corresponding controller $g(\boldsymbol{a})$ .

The non-convexity implies that it is not guaranteed that these iterations will arrive at the global minimum of the bound. Our experience suggests that this is indeed the case and the optimum will typically depend on the initial guess.

It is worth mentioning that, using this bound optimisation procedure, globally stabilising controllers can also be designed. In particular, an SOS globally stabilising controller for an equilibrium point $\boldsymbol{a}_{0}$ is that for which the upper bound optimisation (3.11) has solution $C_{SOS}^{\ast }=\unicode[STIX]{x1D6F7}(\boldsymbol{a}_{0})$ , provided that $\unicode[STIX]{x1D6F7}$ reaches the global minimum on this point. In other words, if the bound can be made tight to $\unicode[STIX]{x1D6F7}(\boldsymbol{a}_{0})$ via control design, then all trajectories of the controlled system must eventually converge to $\boldsymbol{a}_{0}$ , to make the long-time average equal to the bound. Note that here $V$ does not possess the mathematical properties of a Lyapunov function as in the Lyapunov-based method discussed in the previous paragraphs. It is also worth saying that the optimisation stops when the bound cannot be further reduced, resulting in a controller that can still favourably modify the dynamics. Possible reasons for a premature stop include conservativeness of the SOS constraint or a degree of $V$ or $g$ lower than necessary.

4 Application to a fluid flow

In this section the control design methodology described in § 3 is applied to a fluid flow. The mitigation of fully developed vortex shedding, i.e. the nonlinear dynamics of the two-dimensional unsteady wake flow past a circular cylinder at low Reynolds number, $Re=100$ , using a controlled rotary motion of the cylinder, has been selected.

4.1 Numerical set-up

The formulation used to solve the flow problem is based on the Navier–Stokes momentum and continuity equations for a two-dimensional incompressible viscous fluid,

(4.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\frac{1}{Re}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}, & \displaystyle\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$

where $p$ is the reduced pressure and $\boldsymbol{u}=u\boldsymbol{i}+v\boldsymbol{j}$ is the velocity vector defined on a two-dimensional Cartesian space $\boldsymbol{x}=x\boldsymbol{i}+y\boldsymbol{j}$ , centred on the centre of the cylinder, located at $\boldsymbol{x}=(0,0)$ , and oriented such that the $x$ axis is aligned with the free stream, as sketched in figure 2. Normalisation of the governing equations, resulting in (4.1), is done using the cylinder diameter and the free-stream velocity. This yields the standard definition of the Reynolds number as $Re=u_{\infty }{\mathcal{D}}/\unicode[STIX]{x1D708}$ , where ${\mathcal{D}}$ is the cylinder diameter, $u_{\infty }$ is the free-stream velocity and $\unicode[STIX]{x1D708}$ is the kinematic viscosity of the fluid.

Figure 2. Schematic of the problem configuration for the circular cylinder flow. Boundary conditions on the outer domain boundaries are also indicated.

The Navier–Stokes problem (4.1) is solved on a triangular unstructured mesh with a finite volume formulation provided by the open source software OpenFOAM (Jasak, Jemcov & Tukovic 2007). The application icofoam, implementing the well-known PISO algorithm, has been used to solve the velocity–pressure coupling (Ferziger & Perić 2002). Preliminary validation and grid convergence studies, not reported in this paper because of the standard problem type, have been conducted to assess the reliability of the solver, showing good agreement between present and previous numerical results. A mesh of intermediate fineness with size of the elements adjacent to the cylinder equal to 0.02 has been chosen, for a total of approximately 17 000 triangular cells.

The computational domain has the same size as the one used in Bergmann, Cordier & Brancher (2005). It is rectangular and extends for 10 and 20 diameters upstream and downstream of the cylinder, respectively, and spans a total vertical size of 20 diameters in the crossflow. The boundary conditions associated with the problem are also sketched in figure 2. At the inflow, the Dirichlet condition $\boldsymbol{u}=(u_{\infty },0)$ is used for the velocity, while the Neumann condition $\unicode[STIX]{x2202}p/\unicode[STIX]{x2202}x=0$ is used for the pressure. On the upper and bottom boundaries, a free-slip condition is used for the velocity, such that $\unicode[STIX]{x2202}u/\unicode[STIX]{x2202}y=0$ and $v=0$ . A zero-normal-gradient condition is used for the pressure on these two boundaries. On the cylinder surface, the no-slip condition $\boldsymbol{u}=(0,0)$ is enforced, while the standard zero-normal-pressure-gradient condition is used for the pressure. At the outflow boundary, good numerical results, without spurious reflections, were obtained by using a zero-normal-gradient condition for the velocity, i.e. $\unicode[STIX]{x2202}\boldsymbol{u}/\unicode[STIX]{x2202}x=(0,0)$ , while the Dirichlet condition $p=0$ was set to fix uniquely the pressure field.

The time step was constant and equal to $\unicode[STIX]{x0394}t=0.005$ for the mesh used to obtain all the results reported in the rest of the paper. This choice was adopted to achieve satisfactory temporal resolution and a maximum Courant number in the flow field of the order of 0.7.

In the following we will make use of a standard inner product between vector fields, defined as

(4.2) $$\begin{eqnarray}\langle \boldsymbol{v},\boldsymbol{w}\rangle =\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{w}\,\text{d}\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the flow domain, and the associated norm $\Vert \boldsymbol{v}\Vert =\langle \boldsymbol{v},\boldsymbol{v}\rangle ^{1/2}$ .

4.2 Proper orthogonal decomposition and reduced-order modelling

The SOS-based design methodology requires a finite-dimensional representation of the system dynamics, available explicitly as a set of first-order ordinary differential equations with right-hand side being a polynomial function of the state variables and of the control. Spatial discretisation of (4.1), an infinite-dimensional system modelled by partial differential equations, leads formally to such a system, but the extremely large dimension leads to problems that are not tractable numerically, even for moderately complex flows. Specifically, the computational cost of the solution of the SOS problems discussed above, e.g. equation (3.11), grows extremely quickly with the system size, as will be discussed later. Hence, a model reduction strategy is used in this paper to reduce the size of the dynamical system and allow a numerically tractable solution.

We adopt a standard Galerkin projection method, whereby the full dynamics are projected onto a low-dimensional linear subspace, spanned by appropriately selected basis functions. To begin with, the velocity vector field is assumed to be approximated by the ansatz

(4.3) $$\begin{eqnarray}\boldsymbol{u}^{N}(\boldsymbol{x},t)=\overline{\boldsymbol{u}}(\boldsymbol{x})+\unicode[STIX]{x1D6FE}(t)\boldsymbol{u}(\boldsymbol{x})+\mathop{\sum }_{i=1}^{N}a_{i}(t)\boldsymbol{u}_{i}(\boldsymbol{x}).\end{eqnarray}$$

The velocity field is decomposed into a solenoidal steady base flow $\overline{\boldsymbol{u}}(\boldsymbol{x})$ satisfying homogeneous boundary conditions on the cylinder, a ‘control flow’ $\unicode[STIX]{x1D6FE}(t)\boldsymbol{u}(\boldsymbol{x})$ (see e.g. Graham, Peraire & Tang 1999; Kasnakoğlu, Serrani & Efe 2008) used to lift the time-dependent inhomogeneous boundary conditions on the oscillating cylinder surface and to include control via the boundary in the dynamic model, and the weighted sum of $N$ solenoidal vector fields $\boldsymbol{u}_{i}(\boldsymbol{x})$ , the basis functions, which are assumed to form an orthonormal set.

Because the dynamics of a high-dimensional system are compressed into few degrees of freedom, the choice of the basis functions $\boldsymbol{u}_{i}$ is often crucial. Growing interest in model-based control of fluid flow has resulted in different selection strategies, which are far too numerous to discuss here (see e.g. Noack et al. 2003; Barbagallo et al. 2009). In this work we used POD (Sirovich 1987; Berkooz, Holmes & Lumley 1993; Holmes, Lumley & Berkooz 1998) to identify the low-dimensional subspace. The motivating observation for the choice of POD in the current context is that, when data used in the POD algorithm are specifically obtained by sampling the system after the developed regime has established, i.e. any transient has died out, the basis functions describe approximately the axis of inertia of the attractor of the system. Hence, a reduced-order model (ROM) that describes accurately the long-term behaviour of the system is extremely important in the present case, because the focus of the present SOS paradigm in on the estimation and optimisation of a bound for the developed regime.

With the idea of exciting transient flow structures and obtaining a richer snapshot set (Bergmann et al. 2005), a first set of snapshots of the velocity vector field, ${\mathcal{U}}=\{\boldsymbol{u}(\boldsymbol{x},t_{k})\}_{k=1}^{M}$ , is sampled from a DNS in which the angular motion of the cylinder is driven by a random actuation signal. The discrete-time actuation signal is obtained from samples of a zero-mean Gaussian distribution. It is then digitally filtered, such that its power spectrum has zero energy outside the band of reduced frequency $S_{t}=f{\mathcal{D}}/u_{\infty }=[0.1,0.25]$ . The amplitude of the filtered signal is then modulated by a low-frequency mode, $St_{mod}=0.005$ , in order to actuate the flow at different intensities, and it is then normalised to have unitary maximum magnitude, resulting in a standard deviation equal to approximately 0.25. One-third of the total control signal is shown in figure 3. The total duration of the simulation is $T=1000$ , approximately 150 oscillation cycles of the uncontrolled flow, and a total of $M=900$ snapshots is sampled, from $t\geqslant 100$ , at intervals of one non-dimensional time unit. We have verified that such a number of snapshots is sufficient to provide convergence of the second-order statistics associated with POD, such as the energy distribution among individual POD modes.

Figure 3. Time history of the signal used to generate snapshots of the actuated velocity field.

The time-dependent inhomogeneous boundary conditions on the cylinder are lifted from the snapshots by subtracting, with appropriate amplitude, the control function, obtaining the set

(4.4) $$\begin{eqnarray}{\mathcal{U}}^{\prime }=\{\boldsymbol{u}^{h}(\boldsymbol{x},t_{k})=\boldsymbol{u}(\boldsymbol{x},t_{k})-\unicode[STIX]{x1D6FE}(t_{k})\boldsymbol{u}_{c}(\boldsymbol{x})\}_{k=1}^{M}.\end{eqnarray}$$

A radially symmetric solenoidal control function $\boldsymbol{u}_{c}(\boldsymbol{x})$ , with circumferential velocity decaying as $\text{e}^{-\unicode[STIX]{x1D706}(r-0.5)}$ , was used. The decay factor $\unicode[STIX]{x1D706}=8$ was selected from a numerical simulation where the cylinder is oscillated harmonically, in quiescent fluid, with frequency equal to the shedding frequency, and monitoring the decay with $r$ of the amplitude of the velocity fluctuations in the circumferential Stokes flow.

The arithmetic average

(4.5) $$\begin{eqnarray}\overline{\boldsymbol{u}}(\boldsymbol{x})=\frac{1}{M}\mathop{\sum }_{k=1}^{M}\boldsymbol{u}^{h}(\boldsymbol{x},t_{k})\end{eqnarray}$$

is used as the base flow for the ansatz (4.3). Finally, the snapshot set

(4.6) $$\begin{eqnarray}{\mathcal{U}}^{\prime \prime }=\{\boldsymbol{u}^{h}(\boldsymbol{x},t_{k})-\overline{\boldsymbol{u}}(\boldsymbol{x})\}_{k=1}^{M}\end{eqnarray}$$

is used for the POD algorithm. As is common practice, the ‘snapshot’ method of Sirovich (1987) is used.

The selection of the number of basis functions used for the projection, and hence the dimension of the state vector $\boldsymbol{a}$ , is driven by the trade-off between accuracy and cost of computations. The key aspect in this selection is that the computational cost associated with the solution of SOS problems (3.10)–(3.11) increases quite dramatically with the state dimension $N$ and the degrees of the function $V$ and the controller $g$ , $d_{V}$ and $d_{g}$ . For example, the SOS constraint in (3.10) is a polynomial in $N$ variables of degree $d_{V}+1=2d$ , assuming that $\unicode[STIX]{x1D6F7}$ has lower degree than this, and because $\boldsymbol{f}$ is at most quadratic in $\boldsymbol{a}$ for models of incompressible fluid flows. For such a polynomial, the vector of monomials equivalent to $\boldsymbol{v}$ in (2.2) consists of $D=(N+d)!/(N!d!)$ individual terms, whereas the cost for solving the SDP problem in each iteration increases in practice as ${\mathcal{O}}(D^{3})$ (for more details, see Goulart & Chernyshenko 2012, and references therein).

Figure 4. Normalised cumulative energy associated with the POD modes obtained from snapshots sampled from DNS with random actuation. Nine POD modes are selected, capturing 91 % of the total fluctuation kinetic energy associated with the snapshots.

As a compromise between computational cost and model performance, we selected the first nine POD modes for the Galerkin projection, capturing approximately 91 % of the total fluctuation kinetic energy in the snapshots, as illustrated in figure 4, which shows the normalised cumulative energy associated with the POD modes. In addition, this model is augmented with a tenth, shift mode (Noack et al. 2003), a particular mode spanning the direction from the mean flow $\overline{\boldsymbol{u}}(\boldsymbol{x})$ to the unstable, steady and symmetric solution $\boldsymbol{u}_{0}(\boldsymbol{x})$ of (4.1). The symmetric flow is obtained numerically as the steady-state solution using the half upper grid of the original problem, with free-slip boundary condition on the symmetry plane, sufficient to suppress the symmetry-breaking unstable normal mode that grows and saturates into the periodic von Kármán street (Protas & Wesfreid 2002; Bergmann et al. 2005). The shift mode is then constructed as

(4.7) $$\begin{eqnarray}\boldsymbol{u}_{\unicode[STIX]{x1D6E5}}(\boldsymbol{x})=\frac{\overline{\boldsymbol{u}}-\boldsymbol{u}_{0}}{\Vert \overline{\boldsymbol{u}}-\boldsymbol{u}_{0}\Vert },\end{eqnarray}$$

and it is then made orthogonal to the remaining nine POD modes using a Gram–Schmidt procedure.

It is well recognised (Tadmor et al. 2010) that the inclusion of shift modes in Galerkin models of natural and actuated wake flows past a circular cylinder improves transient dynamics over larger changes in base flow. However, a more important result is that inclusion of such a mode results in a dynamical system for which a finite upper bound on the long-time average of energy $\unicode[STIX]{x1D6F7}(\boldsymbol{a})=\boldsymbol{a}^{\text{T}}\boldsymbol{a}/2$ can be found using the procedure presented in § 3.2, similarly to what was observed in Schlegel & Noack (2015). On the other hand, the nine-mode POD-based system does not appear to have such a property, as we have been unable to find an upper bound for the same quantity. Even though the inability to find an upper bound using SOS does not prove that a finite upper bound does not exist, as discussed in § 3.2, it prevents the application of the control methodology proposed in this paper, entirely based on bounds estimation and optimisation.

Standard Galerkin projection is then performed by inserting the expansion (4.3) in (4.1), and setting the inner product with each of the modes to zero in turn. Neglecting the small contribution arising from the projection onto the pressure gradient field, as commonly done for this fluid flow (e.g. Noack et al. 2003; Bergmann et al. 2005), results in the nonlinear system of first-order coupled ordinary differential equations, the ROM:

(4.8) $$\begin{eqnarray}\frac{\text{d}a_{i}}{\text{d}t}=c_{i}+\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D613}_{ij}a_{j}+\mathop{\sum }_{j=1}^{N}\mathop{\sum }_{k=1}^{N}\unicode[STIX]{x1D618}_{ijk}a_{j}a_{k}+m_{i}\frac{\text{d}\unicode[STIX]{x1D6FE}}{\text{d}t}+e_{i}\unicode[STIX]{x1D6FE}+b_{i}\unicode[STIX]{x1D6FE}^{2}+\mathop{\sum }_{j=1}^{N}\unicode[STIX]{x1D60D}_{ij}a_{j}\unicode[STIX]{x1D6FE},\quad i=1,\ldots ,N.\end{eqnarray}$$

The definitions of the coefficients $c_{i},\unicode[STIX]{x1D613}_{ij},\unicode[STIX]{x1D618}_{ijk},m_{i},e_{i},b_{i}$ and $\unicode[STIX]{x1D60D}_{ij}$ arising from the projection are standard and are reported in appendix C. Numerical time integration of the ROM is performed using a standard fourth-order Runge–Kutta scheme with time step equal to  $10^{-3}$ .

4.3 Choice of the cost function

It has been pointed out in the literature (Homescu, Navon & Li 2002) that the choice of the quantity to be minimised by control can sometimes determine the performance of the resulting controller. Hence, several options have been proposed. For instance, in the optimal control approach of Protas & Styczek (2002), using full-order simulations of the Navier–Stokes equations and their adjoint, the chosen cost was the sum of the work needed to resist the drag force and the work needed to control the flow. These two quantities could be computed exactly for that case, but for reduced-order Galerkin-type models, such a level of detail is typically not available, or requires extension of the POD basis to pressure (Bergmann, Bruneau & Iollo 2009). In some works (Graham et al. 1999; Bergmann et al. 2005), the unsteadiness in the wake is typically used as a proxy for drag, and an additional penalisation on the control magnitude is added for regularisation purposes. In the present work, we adopted this formulation, where the cost to be reduced is the domain integral of the kinetic energy of the velocity fluctuations resolved by the ansatz (4.3), plus a penalisation on the control, i.e. the quantity

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}(\boldsymbol{a}(t))={\textstyle \frac{1}{2}}\boldsymbol{a}(t)^{\text{T}}\boldsymbol{a}(t)+R\unicode[STIX]{x1D6FE}^{2}(\boldsymbol{a}(t)),\end{eqnarray}$$

where the orthonormality of the basis functions has been used. The penalisation factor $R$ does not have an immediate physical meaning, but it is used as a design parameter as a means to artificially limit the amplitude of the control. This is necessary because increasingly large control inputs will drive the ROM increasingly far from the region of the phase space where accurate and realistic dynamical behaviour can be expected. As a result, performance in DNS can be affected, as will be shown later.

Figure 5. (a,b) Time histories of states $a_{3}$ and $a_{5}$ from numerical integration of the ROM obtained directly from Galerkin projection (black dashed line) and of the calibrated ROM (red solid line), compared with the time history of the projections of the corresponding POD modes onto the DNS solution. (c) Time histories of system energy for the original and calibrated ROMs, and from projections on the DNS of the uncontrolled flow.

4.4 Model calibration

The 10-mode ROM obtained directly from Galerkin projection is able to represent the dynamics of the full-order system only over a short time scale, i.e. about one shedding cycle, and the long-term behaviour is not correctly represented. This phenomenon is illustrated in figure 5. Figure 5(a,b) shows time histories of the projections of the third and fifth POD modes, respectively, onto the DNS of the uncontrolled flow, i.e. the quantities

(4.10) $$\begin{eqnarray}\tilde{a}_{i}(t)=\langle \boldsymbol{u}_{i}(\boldsymbol{x}),\boldsymbol{u}(\boldsymbol{x},t)-\overline{\boldsymbol{u}}(\boldsymbol{x})\rangle ,\quad i=3,5.\end{eqnarray}$$

These are compared with the time histories of the same quantities obtained from numerical integration of the original ROM, with initial condition $\boldsymbol{a}(0)=\tilde{\boldsymbol{a}}(0)$ (black dashed line). It is clear that the predictions of the model quickly diverge and become essentially useless. The energy of the system $\boldsymbol{a}^{\text{T}}\boldsymbol{a}/2$ (figure 5 c) grows significantly and its long-time average is 12 times larger than the mean resolved energy obtained from projections of the modes onto the DNS solution. The attractor of the ROM is thus significantly different from the projection of the stable limit cycle associated with vortex shedding onto the 10-dimensional phase space. This is a recurrent problem in reducing the order of nonlinear dynamical systems (Cordier et al. 2013) because the high sensitivity to perturbations, such as the truncation of low-energy modes, can have a profound effect after a sufficiently long time (Marion & Temam 1989). A crucial point is that time averages, and bound estimation/optimisation, depend strongly on the geometry of the phase space. It is then very important for an effective application of the proposed control design methodology to have a ROM matching as precisely as possible the long-term behaviour of the original full-order system.

Hence, we apply an eddy-viscosity model calibration scheme, which has become standard practice to correct the effects of unresolved truncated modes (Sirisup & Karniadakis 2004; Couplet, Basdevant & Sagaut 2005; Noack et al. 2008; Bergmann et al. 2009; Protas, Noack & Östh 2015). Following previous work (Cordier et al. 2010), we add to (3.1) a linear calibration term $\unicode[STIX]{x1D613}_{ij}^{c}a_{j}$ , where the matrix $\unicode[STIX]{x1D613}_{ij}^{c}$ has non-zero, initially undetermined, entries only on the main diagonal and the first upper/lower diagonals. Adding the contribution from the two off-diagonals, as opposed to previous work where only the diagonal elements are identified (Galletti et al. 2004), was necessary to achieve satisfactory tracking of the reference limit cycle. Optimal entries are obtained from the solution of the optimisation problem

(4.11) $$\begin{eqnarray}\min _{\unicode[STIX]{x1D613}_{ij}^{c}}~\displaystyle \int _{t_{0}}^{t_{1}}\Vert \boldsymbol{a}(t;\unicode[STIX]{x1D613}_{ij}^{c})-\tilde{\boldsymbol{a}}(t)\Vert ^{2}\,\text{d}t\end{eqnarray}$$

subject to the state equation (4.8) with $\unicode[STIX]{x1D6FE}(t)=0$ , with initial condition $\boldsymbol{a}(t_{0})=\tilde{\boldsymbol{a}}(t_{0})$ . In (4.11), the time integral of the norm of the error between the calibrated ROM trajectory $\boldsymbol{a}(t;\unicode[STIX]{x1D613}_{ij}^{c})$ and the projection of the trajectory of the full-order system onto the selected subspace $\tilde{\boldsymbol{a}}(t)$ is minimised. This trajectory is obtained from numerical simulation of the uncontrolled flow, after transients have died out, in order to force the calibrated ROM to describe correctly the stable limit cycle associated with vortex shedding. A sequence of optimisation problems with increasing $t_{1}-t_{0}$ is formulated to improve convergence of this non-convex problem, with the final interval amounting to approximately 20 shedding cycles. This procedure is not guaranteed to result in successful identification in the general case, but was successful in this case.

Once calibrated, the ROM shows a more realistic behaviour, as the time-averaged energy on the attractor is only 4 % greater than that associated with the limit cycle of the full-order system, as illustrated in figure 5 (solid red lines). However, poor controllability was observed, as opposed to larger models that do not present this behaviour, suggesting that the rotary actuation of the cylinder affects via viscosity the large-scale motions, i.e. the resolved modes, through linear/nonlinear interaction with the truncated modes. To mitigate this poor controllability, we added two additional calibration terms $e_{i}^{c}a_{i}$ and $m_{i}^{c}\,\text{d}\unicode[STIX]{x1D6FE}/\text{d}t$ in (4.8). Optimal values are obtained from an optimisation problem similar to (4.11), where the numerical simulation used to determine the POD modes is used as reference. Although the trajectory of the ROM remains bounded when integrated using the same actuation signal of DNS, it quickly diverges from the reference trajectory and rapidly becomes uncorrelated. Hence, we adopted a more robust multiple-shooting identification scheme (Peifer & Timmer 2007; Protas et al. 2015). The idea is to consider a set of $K$ blocks, each of length $T$ equal to the shedding period, and minimise the sum of all deviations, i.e. solving

(4.12) $$\begin{eqnarray}\min _{e_{i}^{c},m_{i}^{c}}~\mathop{\sum }_{k=1}^{K}\int _{t_{k}}^{t_{k}+T}\Vert \boldsymbol{a}(t;e_{i}^{c},m_{i}^{c})-\tilde{\boldsymbol{a}}(t)\Vert ^{2}\,\text{d}t,\end{eqnarray}$$

subject to the state equation (4.8), with identical notation as in (4.11), and where the initial condition $\boldsymbol{a}(t_{k})=\tilde{\boldsymbol{a}}(t_{k})$ is used for numerical integration of the calibrated ROM over the $k$ th block.

5 Results

The SOS-based methodology discussed above has been used to derive a set of linear state-feedback controllers, i.e. the degree $d_{g}$ has been set to one, with $d_{V}=4$ in all cases, for various penalisation factors. Additional tests have been performed with $d_{V}=6$ with no difference, except for additional computational costs, as all bounds are tight to the actual average from simulation with $d_{V}=4$ . Linear controllers have the form $\unicode[STIX]{x1D6FE}(t)=\sum _{i=1}^{N}k_{i}a_{i}(t)$ , where all the gains have to be identified. The constant term is manually set to zero, i.e. the gains $k_{i}$ are the only decision variables in the SOS calculations. This is done explicitly to avoid naturally occurring spurious control solutions with large non-zero mean rotation, a result probably exploiting unrealistic dynamics described by the ROM far away from the operating regime. It is possible to get rid of the term $m_{i}\,\text{d}\unicode[STIX]{x1D6FE}/\text{d}t$ in (4.8), by noting that $\text{d}\unicode[STIX]{x1D6FE}/\text{d}t=\sum _{i=1}^{N}k_{i}\,\text{d}a_{i}/\text{d}t$ and using the state equation to get

(5.1) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FE}}{\text{d}t}=\frac{1}{1-\displaystyle \mathop{\sum }_{l=1}^{N}k_{l}m_{l}}\mathop{\sum }_{i=1}^{N}k_{i}(c_{i}+\unicode[STIX]{x1D613}_{ij}a_{j}+\unicode[STIX]{x1D618}_{ijk}a_{j}a_{k}+e_{i}\unicode[STIX]{x1D6FE}+b_{i}\unicode[STIX]{x1D6FE}^{2}+\unicode[STIX]{x1D60D}_{ij}a_{j}\unicode[STIX]{x1D6FE}),\end{eqnarray}$$

which is substituted back into the state equation (4.8). This method cannot be used for nonlinear controllers, as the denominator of the fraction in (5.1) would contain an expression in $\boldsymbol{a}$ , making the resulting system non-polynomial in the state variables. A different approach is required as described in appendix D. The degree of the polynomial-type feedback controller can be regarded as a design parameter, and it is just for the sake of simplicity that we consider linear controllers only in this paper, leaving the derivation and testing of nonlinear controllers as future work. It is worth pointing out that, even though the feedback is a linear function of the state, the control design process is aware and exploits the fully nonlinear dynamics of the ROM. No linearisation is performed, unlike in Aleksić-Roessner et al. (2013), who studied a similar feedback control configuration.

Feedback control results are reported for the ROM first, and then for the implementation in DNS.

5.1 Bound estimation and optimisation

An estimate of the long-time-averaged cost (4.9) was obtained by long numerical integrations of the ROM without control, starting from several random initial conditions. All trajectories converged to the same stable limit cycle and the associated long-time-averaged cost was $\overline{\unicode[STIX]{x1D6F7}}^{0}=3.07$ . Trajectories never exhibited blowup, nor converged to a different attractor, although these numerical experiments cannot, of course, be considered as a proof that another stable attractor does not exist in the phase space of the ROM.

The estimation of the upper bound via SOS is performed by trial and error. For a given $C$ we try to find $V$ satisfying the constraint in (3.10). If this is successful in the sense that the resultant SOS decomposition satisfies the feasibility-checking condition (2.4) (Löfberg 2009), we decrease $C$ by $\unicode[STIX]{x1D6FF}C$ , which is 0.01 here, and repeat the trial. In checking the feasibility of the problem, it is important to consider that SOS problems such as (3.10) lead quickly to large SDPs, typically becoming strongly ill-conditioned as the size increases, although the numerical algorithms are based on convex programming. As discussed in detail in Löfberg (2009), the equality constraints associated with (2.2) are only satisfied in the limit of the solution process, as a result of finite-precision arithmetic and various termination criteria.

Figure 6. Performance of linear feedback controllers for various penalisation factors $R$ in closed-loop simulation of the ROM: crosses, upper bound of the long-time-averaged cost; open circles, converged value of the long-time-averaged cost; closed circles, long-time average of the resolved fluctuation kinetic energy. The horizontal line denotes the time average/upper bound for the uncontrolled system.

Table 1. Linear feedback control results for the ROM for different penalisation factors  $R$ .

Several linear controllers have been designed for increasing values of the penalisation factor $R$ . Large values of $R$ have been used, as these lead to better performance in DNS. In figure 6 the performance of these controllers in closed-loop simulation of the ROM is summarised. Long-time averages of the cost are computed from numerical simulations started from an initial condition on the ROM’s limit cycle and by discarding initial transients as control is activated. Figure 6 reports the upper bound $C^{\ast }$ (crosses), the actual time average $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ from simulation (open circles) and the long-time average of the resolved fluctuation kinetic energy $\boldsymbol{a}^{\text{T}}\boldsymbol{a}/2$ (closed circles), with the difference between the latter two quantities being the average cost of control. The horizontal line is the value for the uncontrolled system. Numerical values of the points displayed in figure 6 are also reported in table 1, together with other quantities of interest. The SOS-based control design successfully reduces the upper bound of the system. The reduction is larger for small $R$ , as larger control magnitudes are allowed, as can be deduced by the last column of table 1, which shows the root-mean-square value of the control input. The maximum reduction of the bound is relatively small, i.e. approximately 6 % for $R=50$ ; larger reduction can be found for smaller penalisation factors, although these controllers performed poorly in DNS. A significant part of the total time-averaged cost comes, artificially, from the control. In fact, the time-averaged resolved fluctuation kinetic energy decreases by as much as 13 % for $R=50$ and by 11 % for $R=200$ .

Repeated integration of the controlled ROM, from several random initial conditions, shows that the time average $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ is always below the upper bound, and no instability of the closed-loop system has been observed. The upper bound is tight to the actual average, within the uncertainty of the numerics involved in solution of the SOS problems, as for the estimation of the bound for the uncontrolled case. The upper bound and the average from simulation appear to converge asymptotically to the bound of the uncontrolled system as $R$ increases.

Figure 7. Transient dynamics of the controlled ROM for $R=150$ . (a,c) The trajectory of the ROM projected onto two relevant subspaces. The long-term behaviour of the system is indicated by the uncontrolled and controlled limit cycles. (b) The time histories of the total cost and of the resolved fluctuation kinetic energy. (d) The time history of the control input.

Figure 7 shows the effects of control on the dynamics of the ROM, for the linear controller with $R=150$ , reported as an example, as the other controllers have a qualitatively similar impact on the dynamics. Figures 7(a) and 7(c) show the trajectory of the ROM projected in the $(a_{1},a_{2})$ and $(a_{3},a_{4})$ planes, respectively. The red/blue ‘controlled’/‘uncontrolled’ orbits denote the limit cycle before/after the activation of control. The transient between the two is indicated in light grey. The actuated dynamics converge to a controlled limit cycle, over which the mean resolved kinetic energy is reduced. Under the action of control, the energy of the first two modes decreases, whereas that of modes $a_{3}$ and $a_{4}$ increases slightly. Physically, such a shift is interpreted as a restructuring of the wake, as both pairs of modes correspond spatially to velocity fluctuations oscillating at the shedding frequency. Figure 7(b) shows time histories of the resolved energy (solid black line) and of the total cost (dashed red line). Feedback control is started at $t=10$ . As soon as control is activated, the total cost $\unicode[STIX]{x1D6F7}$ jumps up to approximately 5.5, because the control input, shown in figure 7(c) quickly jumps to approximately 0.15. As anticipated, the penalisation in the control in the cost function does not limit the instantaneous value of the control, as a hard saturation would, but limits only its long-time-averaged contribution.

The resolved kinetic energy decreases substantially in a short transient that takes approximately 10 time units, i.e. just less than two shedding cycles. On the other hand, the total cost takes a longer time to settle to the steady state, approximately 70 time units after activation of control, because the peak-to-peak variation of the control input $\unicode[STIX]{x1D6FE}(t)$ decreases slowly during the transient.

5.2 Feedback control in DNS

The four linear controllers derived are implemented in DNS. Because the governing equations are marched forward in time, at the beginning of the $k$ th time step, at time $t_{k}$ , the current state vector is obtained from the projections of the POD modes on the current fluctuating velocity field as

(5.2) $$\begin{eqnarray}a_{i}(t_{k})=\langle \boldsymbol{u}_{i}(\boldsymbol{x}),\boldsymbol{u}(\boldsymbol{x},t_{k})-\overline{\boldsymbol{u}}(\boldsymbol{x})-\unicode[STIX]{x1D6FE}(t_{k})\boldsymbol{u}_{c}(\boldsymbol{x})\rangle ,\quad i=1,\ldots ,N.\end{eqnarray}$$

The control action $\unicode[STIX]{x1D6FE}(t_{k})$ , the tangential velocity on the cylinder surface, is calculated from the control law (2.2) and is then set as constant boundary condition for the time step $t_{k+1}-t_{k}$ .

Figure 8 shows time histories of three key quantities obtained from DNS of the closed-loop system: the fluctuation kinetic energy resolved by the Galerkin expansion (4.3),

(5.3) $$\begin{eqnarray}K_{\boldsymbol{u}_{G}^{\prime }}(t)=\frac{1}{2}\Vert \boldsymbol{u}_{G}^{\prime }(\boldsymbol{x},t)\Vert ^{2}=\frac{1}{2}\left\Vert \mathop{\sum }_{i=1}^{N}\tilde{a}_{i}(t)\boldsymbol{u}_{i}(\boldsymbol{x})\right\Vert ^{2}=\frac{1}{2}\tilde{\boldsymbol{a}}(t)^{\text{T}}\tilde{\boldsymbol{a}}(t),\end{eqnarray}$$

obtained from the projections (5.2), in figure 8(a); the total fluctuation kinetic energy,

(5.4) $$\begin{eqnarray}K_{\boldsymbol{u}^{\prime }}(t)={\textstyle \frac{1}{2}}\Vert \boldsymbol{u}(\boldsymbol{x},t)-\overline{\boldsymbol{u}}(\boldsymbol{x})-\unicode[STIX]{x1D6FE}(t)\boldsymbol{u}_{c}(\boldsymbol{x})\Vert ^{2},\end{eqnarray}$$

in figure 8(b); and the kinetic energy of the residual fluctuation $\boldsymbol{u}_{r}(\boldsymbol{x},t)=\boldsymbol{u}^{\prime }(\boldsymbol{x},t)-\boldsymbol{u}_{G}^{\prime }(\boldsymbol{x},t)$ ,

(5.5) $$\begin{eqnarray}K_{\boldsymbol{u}_{r}^{\prime }}(t)=K_{\boldsymbol{u}^{\prime }}(t)-K_{\boldsymbol{u}_{G}^{\prime }}(t),\end{eqnarray}$$

normalised with the total fluctuation kinetic energy $K_{\boldsymbol{u}^{\prime }}$ , in figure 8(c). Note that the cost of the control $R\unicode[STIX]{x1D6FE}(t)^{2}$ is not added to panel $(a)$ . The lower panels (figure 8 d–f) show the same quantities in the interval $t\in [110,155]$ , in the initial transient after activation of feedback control at $t=112$ .

Figure 8. Performance of linear feedback controllers in DNS: $(a)$  time history of fluctuation kinetic energy resolved by the ansatz (4.3);  $(b)$  time history of the total fluctuation kinetic energy; $(c)$  time history of the unresolved residual energy, normalised with the total fluctuation kinetic energy. (d–f) The same quantities as (a–c), respectively, in the interval $t\in [110,155]$ .

As soon as control is activated, the resolved and total kinetic energy decrease substantially, in a transient lasting for approximately 10–12 time units, similarly to that exhibited by the ROM in figure 7. The initial time rate of change of the energy and the maximum reduction are larger for smaller penalisation factors, as the control is more aggressive. Subsequently, the cost remains approximately constant for a short period, in the interval $t\in [125,140]$ . For $R=200$ , $K_{\boldsymbol{u}_{G}^{\prime }}$ has an average value in this window roughly equal to that obtained in simulation of the ROM in figure 6. For the lower penalisation factors, the reduction of the resolved energy is larger than that obtained in simulation of the ROM, suggesting that the ROM significantly underestimates the effects of control on the dynamics. The reduction of the resolved and, most importantly, of the total fluctuation energy suggests that control design has successfully identified the control mechanism to attenuate vortex shedding.

In figure 8 $(f)$ , the fraction of unresolved kinetic energy grows significantly in this first interval, from a value of approximately 4 %, up to approximately 20 % for the smaller $R$ . This shows that, under the effects of control, the full-order system explores regions of the high-dimensional phase space not included in the initial low-dimensional subspace chosen for the projection, especially for larger control inputs. Taking into account the slow deformation of the wake structure unaccounted for in the original POD basis, using, for example, deformable modes (Tadmor et al. 2011) or updating the mode set (see Bergmann et al. 2009, and references therein), might be beneficial to limit this behaviour and achieve improved performance.

After these initial stages, the character of the solution depends strongly on the penalisation $R$ . For $R=150$ and $200$ , the long-term behaviour of the system is a controlled limit cycle with a reduced fluctuation kinetic energy, as predicted by the ROM in figure 7. By contrast, for the two lower penalisation factors, the structure of the long-term behaviour is significantly different. The time history of the total fluctuation kinetic energy (figure 8 $b$ ), undergoes a periodic low-frequency, large-amplitude modulation, with a period of approximately 18 time units, not clearly visible from the resolved energy in figure 8 $(a)$ .

Figure 9. The top six panels show snapshots of vorticity from DNS of the controlled flow with linear controller with $R=50$ . The bottom panel shows the time history of the total fluctuation kinetic energy. Vertical lines denote the time instants at which snapshots are extracted, at $t=112,121,136,142,156$ and $162$ .

Insight into this phenomenon can be gained by analysing in more detail the behaviour of the system from approximately $t=135$ onwards, indicated with a dashed vertical segment in figure 8( $d$ – $f$ ). The total kinetic energy, panel $(e)$ , and the normalised unresolved energy, panel $(f)$ , increase significantly for $R=50$ and 100. The resolved fluctuation energy is practically constant in the interval $t\in [135,140]$ , panel $(d)$ , and only when this growth saturates does the quantity $K_{\boldsymbol{u}_{G}^{\prime }}$ begin to increase. The physical mechanism responsible for this behaviour is illustrated in figure 9, for $R=50$ . The upper panels of figure 9 show six snapshots of the vorticity field, with the colour map clipped at $\pm 3$ to visualise the structure of the actuated wake, although the maximum vorticity magnitude can be as high as 25 in the boundary layer on the cylinder. The bottom panel shows the time history of the total fluctuating kinetic energy, also reported in figure 8 $(b)$ . The vertical lines denote the times $t_{i}$ at which the snapshots are extracted.

In the initial transient after activation of control, between $t=112$ and $t=121$ (snapshots ( $a$ ) and ( $b$ ) in figure 9), the controlled rotation of the cylinder reorganises the generation and dynamics of vorticity in the near wake. The roll-up of the two shear layers is delayed and the fluctuation energy decreases steadily. Shortly after time $t_{b}$ , the fluctuation energy stops decreasing and during the interval $[t_{b},t_{c}]$ the wake locks onto an actuated limit cycle, with reduced $K_{\boldsymbol{u}^{\prime }}$ . The structure of the wake in this regime, panel $(c)$ , is significantly different from the unactuated wake of panel $(a)$ . It is narrower, especially in $4\lesssim x\lesssim 9$ , and the streamwise separation between the structures is shorter. Shortly after $t_{c}=136$ , the fluctuation kinetic energy grows rapidly. This event is connected to the breakdown of the wake structure of panel $(c)$ in figure 9, arising as a large-scale flapping of the entire near-wake flow. This effect might be connected to the growth of an instability of this wake structure, or it might be simply induced by the feedback. After $K_{\boldsymbol{u}^{\prime }}$ peaks, the restructuring of the wake into the original uncontrolled state enables the control to reduce the fluctuation energy, panel $(e)$ , although the same break-up as observed in snapshot $(c)$ occurs shortly after $t_{e}$ . This mechanism repeats indefinitely originating the low-frequency modulation visible in figure 8.

Although the total fluctuation kinetic energy is successfully reduced in DNS of the closed-loop system, the long-time average of the total cost $\unicode[STIX]{x1D6F7}(\tilde{\boldsymbol{a}}(t),\tilde{\unicode[STIX]{x1D6FE}}(t))$ , as defined in (4.9), is not decreased. This result is shown in figure 10 $(a)$ , whereas panel $(b)$ contains a zoom of the same quantity at the initial stages of the simulation, when control is activated. The total cost initially spikes at quite large values, as the control input is quite intense, similarly to what is observed for the ROM in figure 7 for $R=150$ . As control modifies the wake structure, the fluctuations of the cost decrease substantially, below the reference value of the uncontrolled system (horizontal line) approximately in the interval $t\in [120,125]$ . However, the loss of control performance described above in figure 9 results in a strong increase of the total instantaneous cost, especially for the two lower penalisation factors. As a result, the long-time-averaged cost $\overline{\unicode[STIX]{x1D6F7}}^{\ast }$ , reported in table 2, is above the reference value of the uncontrolled system, $\overline{\unicode[STIX]{x1D6F7}}^{0}=2.95$ , also for the two larger penalisation factors.

Figure 10. Time history of the total cost $\unicode[STIX]{x1D6F7}$ , as defined in (4.9). Panel $(b)$ shows a detail of the same time trace, in a small time interval at the early stages of the simulation, indicated by the rectangle in panel  $(a)$ .

The ROM results, table 1, show that, for $R=200$ , the percentage reduction of the cost is quite small, approximately 0.5 %, as a rather large contribution comes from the control cost. In DNS, the same controller results in an increase of the total cost of approximately 5 %. The discrepancy between these two values is certainly within the accuracy of the ROM in describing the effects of actuation on the full-order dynamics.

A physically meaningful quantity is the long-time average of the total power spent to sustain the motion of the cylinder (Bergmann, Cordier & Brancher 2006). This quantity, expressed per unit of span, is the sum of the power $P_{D}$ spent to move the cylinder at speed $u_{\infty }$ against the drag force $D$ and the power required to control the flow, i.e. the power $P_{M}$ required to rotate the cylinder at angular speed $\dot{\unicode[STIX]{x1D703}}$ (positive when anticlockwise) against the viscous torque $M$ exerted by the fluid on the cylinder (positive when it induces a clockwise rotation).

Table 2. Linear feedback control results in DNS for different penalisation factors $R$ . For the uncontrolled system, $\overline{\unicode[STIX]{x1D6F7}}^{\ast }=\overline{\unicode[STIX]{x1D6F7}}^{0}=\overline{\boldsymbol{a}^{\text{T}}\boldsymbol{a}}/2$ . PSR $=$ power saving ratio.

The drag and viscous torque are determined by the dimensional pressure and viscous stress surface distributions $p(\unicode[STIX]{x1D703})$ and $\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D703})$ as

(5.6a,b ) $$\begin{eqnarray}D=-\int _{0}^{2\unicode[STIX]{x03C0}}p(\unicode[STIX]{x1D703})\cos (\unicode[STIX]{x1D703})+\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D703})\sin (\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703},\quad M=-\frac{{\mathcal{D}}}{2}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$

where the viscous stress on the surface arises from the distribution of the tangential velocity $u_{\unicode[STIX]{x1D703}}$ as

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D703})=\unicode[STIX]{x1D707}\left.\!\frac{\unicode[STIX]{x2202}u_{\unicode[STIX]{x1D703}}}{\unicode[STIX]{x2202}r}\right|_{{\mathcal{D}}/2},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is the dynamic viscosity.

The first contribution to the total power spent is then simply

(5.8) $$\begin{eqnarray}P_{D}=Du_{\infty }={\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}u_{\infty }^{3}{\mathcal{D}}C_{D},\end{eqnarray}$$

whereas the second reads as

(5.9) $$\begin{eqnarray}P_{M}=M\dot{\unicode[STIX]{x1D703}}=2M\unicode[STIX]{x1D6FE}u_{\infty }{\mathcal{D}}=\unicode[STIX]{x1D70C}u_{\infty }^{3}{\mathcal{D}}C_{M}\unicode[STIX]{x1D6FE},\end{eqnarray}$$

in which $C_{D}$ and $C_{M}$ are the coefficients of drag and moment, and $\unicode[STIX]{x1D6FE}$ is the normalised surface velocity as introduced above. Note that $P_{M}$ is positive when the cylinder transfers energy to the fluid and negative otherwise. In non-dimensional terms, the total power spent is then expressed by the total power coefficient

(5.10) $$\begin{eqnarray}C_{P}=\frac{P_{D}+P_{M}}{\unicode[STIX]{x1D70C}u_{\infty }^{3}{\mathcal{D}}/2}=C_{D}+2C_{M}\unicode[STIX]{x1D6FE}.\end{eqnarray}$$

The time histories of the total power coefficient obtained from DNS of the closed-loop system are reported in figure 11, for the four controllers derived, with red dashed lines. The drag coefficient is also reported as a black solid line, for reference. The difference between the two, i.e. $C_{P}-C_{D}$ , is the normalised energy per unit time and unit span required to actively control the flow.

Figure 11. Time histories of the total power and drag coefficients, $C_{P}(t)$ and $C_{D}(t)$ , obtained from closed-loop DNS for the four controllers, with $R=50,100,150$ and 200, in panels $(a)$ , $(b)$ , $(c)$ and $(d)$ , respectively. The difference between the two is the energy per unit time transferred to the fluid by the control, i.e. the power spent for actuation.

For $R=50$ and $R=100$ , figure 11( $a$ , $b$ ), the drag exhibits a low-frequency modulation similar to the total fluctuation kinetic energy, with the ‘valleys’ of these two quantities matching fairly well. The drag minima can be as low as 1.28, suggesting that the control design methodology is indeed effective, although performance is periodically lost, as discussed above. Interestingly, the drag coefficient associated with the steady laminar solution $\boldsymbol{u}_{0}$ is, in our set-up, 1.14. As a result, the drag coefficient reduction, compared to the drag associated with vortex shedding (Protas & Styczek 2002), can be, instantaneously, as large as 46 %. The time-averaged drag is of course higher, as reported in table 2. A maximum percentage drag reduction, normalised with the drag coefficient of the uncontrolled flow, of 4.6 % has been obtained, at $R=100$ . This value is not as high as in previous closed-loop control studies on this same configuration. Using optimal control theory, Protas & Styczek (2002) and more recently Flinois & Colonius (2015) have achieved drag reductions of 7 % at $Re=75$ and 15 % at $Re=150$ , and 19 % at $Re=100$ , respectively. However, in these two works, Navier–Stokes equations were used directly for control design, and not a reduction thereof, enabling an effective control strategy to be found. We believe that developing controllers on larger and more accurate ROMs, that correctly describe the change in dynamics as control is activated, will result in increased performance.

On some occasions, the total power coefficient is lower than the drag itself, because the product $C_{M}(t)\unicode[STIX]{x1D6FE}(t)$ is negative. These events indicate that the cylinder is being driven by the viscous torque, corresponding to a passive mechanism where the flow exerts a net work on the cylinder. Nevertheless, this product is positive for most of the time, and indicates that the control is actively manipulating the flow.

For the two larger penalisation factors, the long-term cost of the control is extremely small, with peaks of $C_{P}-C_{D}$ not exceeding 0.002, practically invisible in figure 11( $c$ , $d$ ). The control strategy identified is quite efficient, because, in the long term, a small amount of power is actively spent to reduce the total power by a significantly larger amount. Following Protas & Styczek (2002) and Bergmann et al. (2005), and based on the definition (5.10), the efficiency can be quantified by the power saving ratio (PSR)

(5.11) $$\begin{eqnarray}\text{PSR}=\frac{\overline{C_{P}}^{u}-\overline{C_{P}}^{c}}{2\overline{C_{M}\unicode[STIX]{x1D6FE}}},\end{eqnarray}$$

i.e. the ratio between the power saved and the power spent for control, in a time-averaged sense, where the superscripts $u$ and $c$ indicate the uncontrolled and controlled cases, respectively. The PSR, reported in table 1, is remarkably large for the two higher penalisation factors, when the feedback control operates close to the design point for which it was constructed. For comparison, Protas & Styczek (2002) obtained a PSR equal to 51 at $Re=150$ , and 122 at $Re=75$ , using optimal control theory in a predictive setting. Various other open-loop-type control approaches, where the cylinder is oscillated harmonically at an optimal frequency and amplitude, are significantly less efficient (Bergmann et al. 2005). This is due to the fact that in the present case the feedback controller trades inexpensive control power, scaling directly with the square of the control amplitude and inversely with the square root of the Reynolds number (Bergmann et al. 2006), with precious propulsion power mainly associated with the pressure drag.

6 Discussion and conclusions

The main contribution of the present paper is the development of a novel feedback control design paradigm for incompressible fluid flows. It applies to finite-dimensional dynamical systems given as a set of first-order nonlinear ordinary differential equations, with the right-hand side being a polynomial function in the state variables and in the controls. Galerkin-type models of incompressible fluid flows, obtained from projection of the governing equations on a finite low-dimensional subspace, have exactly this form.

This paradigm of control is rooted on recent advances in control theory and optimisation over polynomials, commonly known as SOS methods. At the core, these methods leverage computationally efficient approaches to construct positive polynomial functions, by formulating and solving convex SDPs.

The key distinguishing features are that (i) the long-term behaviour of the system, the permanent regime, is central in the design stage, i.e. long-time averages of fluctuating quantities can be optimised by control design, and that (ii) the nonlinearity is taken directly into account in the design process. Furthermore, the present SOS-based scheme allows the design of polynomial-type feedback controllers of arbitrary degree, hence is not limited to the linear case discussed here. Further research is required to understand if nonlinear feedback control can considerably improve performance.

We have numerically investigated the problem of mitigating the kinetic energy of velocity fluctuations in the unsteady wake of a circular cylinder at $Re=100$ , in the laminar regime, via controlled rotary motions of the surface, in a full-information state-feedback arrangement. A 10-mode POD–Galerkin ROM of the actuated wake flow was derived. A crucial element is that the phase space of the ROM should host an attractor whose structure is as similar as possible to that of the full-order system when projected on the low-dimensional subspace. This necessity arises from the fact that bound estimation and optimisation via control design target explicitly the attractor of the system and control performance probably increases if the long-term behaviour of the reduced- and full-order systems are similar. From this perspective, model calibration schemes that ensure long-term stability of the reduced system, and similarity of attractors, at least in a statistical sense, are desirable (see e.g. Östh et al. 2014, and references therein).

Linear state-feedback controllers were derived using the ROM, using a penalisation on the control as a design parameter, and were implemented in DNS. These controllers effectively decreased the ‘size’ of the limit cycle associated with vortex shedding and significantly reduced the long-time average of the total fluctuation kinetic energy, as well as the time-averaged drag coefficient. The feedback system was energetically efficient, as the power saved per unit control power spent was in the range of 75–80. For lower values of the penalisation factor, the greater control input resulted in better performance just after activation of control, but eventually performance worsened significantly. This is not a limitation of the present SOS-based scheme, but rather is driven by the POD–Galerkin modelling strategy used, which is known to lack robustness. We expect that improvements in the modelling strategy will result in increased performance in DNS.

Currently, the main drawback of this methodology is the unfavourable scaling of the computational cost with the size of the system $N$ , and the degree of the polynomial function $V$ . Limitations of existing computational tools to preprocess the polynomial inequalities and solve the associated SDPs currently limit the methodology to dynamical systems of size not greater than approximately 10–20, and the size reduces considerably if high-degree polynomials $V$ are used. However, these methods are rather novel and the available computational software tools are designed for generality. Hence, a large number of optimisations can be introduced by specialising these tools to the peculiarities of hydrodynamic-type systems. A few illustrative instances, towards which future efforts will be devoted, are the exploitation of more efficient SDP solvers, sparsity patterns on the right-hand side of the dynamical system as well as in the structure of the tunable function  $V$ .

It is worth pointing out that the improvement in our ability to solve the relevant SDP problems necessary for achieving better results might actually be far less than it might seem. In the case of global stability analysis, the scalability issue was successfully dealt with in Huang et al. (2015) using the uncertain system method proposed in Goulart & Chernyshenko (2012). This approach can be extended to certain problems of flow control. It allows one to construct storage functionals for averaged parameters of systems governed by partial differential equations, that is, systems with an infinite number of degrees of freedom, while solving the SDP problems corresponding to only a limited number of degrees of freedom. The method can also be used to reduce the effective number of degrees of freedom in a finite-dimensional dynamical system. Further details can be found in the cited papers. Here, we only note that, when the bounds are constructed for an uncertain system, the corresponding SDP problems have to deal with twice as many independent variables as in the case of a standard (certain in our terminology) system. Hence, the way forward is to construct a better ROM with, for example, 40 modes, reduce it to an uncertain system with, for example, 15 modes, and then build a controller solving an SDP corresponding to 30 independent variables. The required increase in the quality of the SDP solver then corresponds to only three times increase in the number of independent variables as compared to the case considered in the present paper. This might indeed become possible in the foreseeable future.

A second limitation of the methodology is that the controller is formally guaranteed to reduce only the upper bound of a long-time-averaged cost function. In a particular realisation of the controlled flow, the actual time average might not decrease. The essential motivation, discussed at length in this paper, is that control design targets the attractor for which the time-averaged cost is the largest, i.e. the one with which the upper bound is associated. If the system has, for instance, two different attractors, with separate basins of attraction, the control of a trajectory started from an initial condition not in the basin of the attractor associated with the bound might result in worsened performance. From this perspective, the controller is guaranteed to reduce the time average only in the worst-case scenario, which might be different from the most likely scenario. In practice, this limitation might be less important than it appears here, although it represents a possibility in the general case.

A further observation is that in this paper we have investigated the control of a dynamical system for which the permanent regime is given by a trivial attractor, i.e. a stable periodic orbit. However, real turbulent flows usually have extremely complicated attractors, and the long-term behaviour is usually chaotic. Reduced-order modelling and long-time-average cost control of such a type of system via SOS optimisation would be more interesting, but also would induce more difficulties. We would like to address them in our future work.