Enrico Deusebio
email: deusebio@mech.kth.se

Pierre Augier

E Lindborg

Third-order structure functions in rotating and stratified turbulence: a comparison between numerical, analytical and observational results

Keywords: atmospheric flows, turbulent flows, turbulence theory

HAL is

Introduction

Statistical analysis of atmospheric measurements show that horizontal spectra of both wind and temperature in the upper troposphere and lower stratosphere follow a k -5/3 power law over more than two decades, at wavenumbers corresponding to wavelengths between 2 and 500 km [START_REF] Nastrom | A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft[END_REF]. Such a power law is a strong indication that there is an energy cascade, that is a systematic transfer of energy between different scales by nonlinear interactions. Either there is positive spectral kinetic energy (KE) flux, Π K , which means that energy flows from small to large wavenumbers, as in the downscale energy cascade of three-dimensional isotropic turbulence (Kolmogorov 1941b), or there is a negative spectral energy flux, which means that energy flows in the opposite direction, as in the upscale energy cascade of two-dimensional turbulence [START_REF] Kraichnan | Inertial ranges in two-dimensional turbulence[END_REF] or quasi-geostrophic (QG) turbulence [START_REF] Charney | Geostrophic turbulence[END_REF]. A measurement of the energy spectrum is not sufficient to discriminate between these two types of cascades since they both exhibit an energy spectrum of the form E(k) ∼ |Π K | 2/3 k -5/3 . [START_REF] Dewan | Stratospheric wave spectra resembling turbulence[END_REF] hypothesized that the dynamics producing the mesoscale k -5/3 power law in the upper troposphere and lower stratosphere is a downscale energy cascade of interacting gravity waves, while [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF] hypothesized that the dynamics is an upscale energy cascade of layerwise two-dimensional turbulence. [START_REF] Lindborg | Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence?[END_REF] suggested that measurements of third-order velocity structure functions could give an answer to the question regarding the direction of the energy cascade. To obtain such an answer one should use a relation which is of similar form to the so-called 'four-fifths law' of isotropic three-dimensional turbulence (Kolmogorov 1941a)

δv L δv L δv L = -4 5 Π K r, (1.1)
where • denotes averaging and δv L (r) = v Lv L is the difference between the longitudinal velocity component at two points separated by a vector r, where the longitudinal component is the projection onto the unit vector n = r/r. The mean energy dissipation per unit mass, ε, may be substituted for Π K in (1.1), since they are equal when the cascade is downscale. The four-fifths law is the three-dimensional incompressible and isotropic solution to a more general flux relation which can be written as

∇ • J K = -4Π K ,
(1.2)

where J K = δvδv • δv . The two-dimensional isotropic solution to this relation can be written as [START_REF] Lindborg | Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence?[END_REF])

δv L δv L δv L + δv L δv T δv T = -2Π K r, (1.3)
where δv T is the transverse velocity increment, that is an increment of a velocity component which is perpendicular to n. By incompressibility one can also derive the relation δv L δv L δv L = -3 2 Π K r, (1.4) for two-dimensional turbulence. Comparing the three-dimensional relation (1.1) with the two-dimensional relations (1.3) and (1.4) we see that they have the same form.

The third-order structure function is proportional to the spectral energy flux and is linear in the separation distance. However, there is a very important difference. In three-dimensional turbulence the spectral energy flux is positive while it is negative in two-dimensional turbulence. Therefore, an upscale energy cascade in the mesoscales would mean that the third-order structure function should display a positive linear dependence on r. On the other hand, a negative linear dependence on r would be an indication of a downscale energy cascade. Cho & Lindborg (2001) measured the third-order structure function in the lower stratosphere, with r lying in the horizontal

∇ • J K = -4Π K + 2 δbδw , (1.5) ∇ • J P = -4Π P -2 δbδw .
(1.6)

Here, b = θ g/θ 0 is the buoyancy, θ and θ 0 are the fluctuating and mean potential temperature respectively, g is the acceleration due to gravity, w is the vertical velocity, Π P is the spectral flux of APE and J P = δvδbδb /N 2 , where N is the Brunt-Väisälä frequency. The last term in (1.5), which appears with opposite sign in (1.6), represents energy exchange between KE and APE. In a range of scales where there is a selfsimilar energy cascade of universal character the ratio between KE and APE should be scale independent and this term should therefore be negligible. Neglecting this term and introducing the assumption of axisymmetry the relations (1.5) and (1.6) can be written as

1 r h ∂ ∂r h (r h δu L δv • δv ) + ∂ ∂r z δwδv • δv = -4Π K , (1.7) 1 r h ∂ ∂r h (r h δu L δbδb /N 2 ) + ∂ ∂r z δwδbδb /N 2 = -4Π P , (1.8)
where r h is the radial coordinate, that is the projection of r on the horizontal plane, r z is the axial coordinate and the velocity is decomposed as v = u + we, with u = u L n + u T t the horizontal velocity. Here, t and e are the azimuthal and axial unit vectors, respectively. Taking r to be in the horizontal plane and assuming that all terms including the vertical velocity can be neglected these relations can be integrated to give

δu L δu L δu L + δu L δu T δu T = -2Π K r, (1.9) δu L δbδb /N 2 = -2Π P r, (1.10)
where we have changed the name of the radial coordinate from r h to r, which is natural when r lies in the horizontal plane. The relation (1.9) is identical to the twodimensional relation (1.3) which was used by Cho & Lindborg (2001). [START_REF] Lindborg | Determining the cascade of passive scalar variance in the lower stratosphere[END_REF] also calculated the third-order velocity-temperaturetemperature structure function, δu L δTδT . In figure 2 we show this function divided by 2r, plotted versus separation distance. In the interval 20 km < r < 400 km there is a range where δu L δTδT /(2r) is negative and constant. In fact, this result is much cleaner than the corresponding result for the velocity structure function seen in figure 1, and the range with a negative linear dependence is much wider. We will now interpret this result in the light of relation (1.10), an interpretation that was not given by [START_REF] Lindborg | Determining the cascade of passive scalar variance in the lower stratosphere[END_REF]. The aircraft data are taken at constant-pressure surfaces. Therefore, changes in the potential temperature θ are only due to changes in absolute temperature T. On a such surface we therefore have

δu L δbδb = δu L δTδT g 2 T 2 0 , (1.11)
where T 0 is the mean temperature. Using relations (1.10) and (1.11), the typical lower stratospheric values N 2 = 4 × 10 4 s -2 and T 0 = 220 K [START_REF] Charney | Propagation of planetary-scale disturbances from the lower into the upper atmosphere[END_REF], and the constant value 3.9 × 10 -6 K 2 s -1 from figure 2, we can estimate the spectral flux of APE as Π P ≈ 2 × 10 -5 m 2 s -3 .

(1.12)

With the estimated value of Π K we would have Π K /Π P ≈ 3. The interpretation of the negative linear range of the third-order structure functions from the lower stratosphere, seen in figures 1 and 2, as evidence of a downscale energy cascade has gained considerable support from theoretical and numerical advances in strongly stably stratified dynamics, although recent experiments on electromagnetically forced turbulence in rotating thick layers have also revived the hypothesis of an upscale energy transfer [START_REF] Xia | Upscale energy transfer in thick turbulent fluid layers[END_REF]. Both scaling analysis [START_REF] Billant | Self-similarity of strongly stratified inviscid flows[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF], idealized numerical simulations [START_REF] Laval | Forced stratified turbulence: Successive transitions with Reynolds number[END_REF][START_REF] Riley | Dynamics of turbulence strongly influenced by buoyancy[END_REF][START_REF] Waite | Stratified turbulence dominated by vortical motion[END_REF][START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF][START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]) and experiments [START_REF] Augier | Stratified turbulence forced with columnar dipoles[END_REF][START_REF] Augier | Experimental study on stratified turbulence forced with columnar dipoles[END_REF] show that there exists a strongly stratified turbulent regime at horizontal scales which are considerably larger than the Ozmidov length scale, l 0 = ε 1/2 /N 3/2 . The conditions for such a regime to exist are that the horizontal Froude number F h = U/(LN), where U is a characteristic horizontal velocity and L a characteristic horizontal length scale, is small and that the buoyancy Reynolds number Re = ε/(νN 2 ), where ν is the kinematic viscosity, is large. Due to the strong intermittency of stratified turbulence, there is a strong variation of the Ozmidov length scale in the lower stratosphere. A typical value can be estimated to be around 10 m [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]. Classical isotropic turbulence, where horizontal and vertical velocity fluctuations are of the same order of magnitude, can be observed at scales below the Ozmidov length scale. At larger scales, vertical velocity fluctuations are strongly damped by buoyancy forces.

Yet, in this regime there exists a turbulent energy cascade which is not strictly twodimensional, as initially suggested by [START_REF] Lilly | Stratified turbulence and the mesoscale variability of the atmosphere[END_REF], since the nonlinear advection terms in the momentum equation including the vertical velocity are of the same order of magnitude as the corresponding horizontal advection terms, that is

|w∂ z u| ∼ |u∂ x u|, (1.13)
where x is a horizontal coordinate and u the corresponding horizontal velocity component. The reason why this can be true is that the vertical length scale of the turbulent motions decreases with increased degree of stratification, so that |∂ z u| increases at the same rate as |w| decreases when the stratification becomes strong. In this regime, there is a downscale energy cascade of KE, which is totally dominated by the horizontal velocity component, as well as APE. The estimate Π K /Π P ≈ 3 which we have made on the basis of (1.9), (1.10) and the stratospheric data presented in figures 1 and 2 is consistent with results from numerical simulations of stratified turbulence [START_REF] Lindborg | The energy cascade in a strongly stratified fluid[END_REF]). It has also been demonstrated that the strongly stratified turbulence regime can be obtained in the presence of system rotation [START_REF] Lindborg | The effect of rotation on the mesoscale energy cascade in the free atmosphere[END_REF][START_REF] Waite | The transition from geostrophic to stratified turbulence[END_REF], provided that the rotation is not too strong. Recently, [START_REF] Vallgren | A possible explanation of the atmospheric kinetic and potential energy spectra[END_REF] and [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF] carried out numerical simulations of strongly stratified turbulence at strong system rotation, that is at small Rossby number Ro = U/f L, where f is the Coriolis parameter. They showed that even for very small Ro a substantial part of the injected energy cascades towards small scales. The amount of energy going downscale increases with Ro, i.e. as the rotation rate decreases. This leakage of energy from the large-scale QG flow leads at smaller scales to a strongly stratified turbulent regime with a downscale cascade of KE and APE, displaying k -5/3 energy spectra. At large scales, however, there is a QG regime displaying a k -3 energy spectrum. The superposition of the k -5/3 energy spectrum of stratified turbulence with the k -3 spectrum of QG turbulence gives a compound spectrum which is similar to that measured by [START_REF] Nastrom | A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft[END_REF] in the atmosphere.

An issue which is still debated is what role gravity waves have in the downscale energy cascade, producing the Nastrom-Gage spectrum. As stated previously, it was hypothesized by [START_REF] Dewan | Stratospheric wave spectra resembling turbulence[END_REF] that the downscale cascade is a process in which gravity waves are successively transferring their energy to modes with shorter wavelengths, by nonlinear interaction. [START_REF] Bartello | Geostrophic adjustment and inverse cascades in rotating stratified turbulence[END_REF] argued that the most important triad interaction transferring energy downscale is the type where two wave modes interact with a vortical mode, with the last one acting as catalyst. This picture has recently been advanced by [START_REF] Waite | The mesoscale kinetic energy spectrum of a baroclinic life cycle[END_REF] and [START_REF] Waite | Mesoscale energy spectra of moist baroclinic waves[END_REF] who carried out simulations of an idealized baroclinic wave cycle, showing the spontaneous development of gravity waves which caused the shallowing of the spectrum. On the other hand, it has also been argued that waves play a minor role in the downscale energy cascade. [START_REF] Lindborg | Horizontal wavenumber spectra of vertical vorticity and horizontal divergence in the upper troposphere and lower stratosphere[END_REF] used structure functions calculated from aircraft data to make a rotational-divergent decomposition of the energy spectrum and found that rotational energy is somewhat larger than divergent energy, from which it can be argued that gravity waves cannot be the only source of energy in the k -5/3 -spectrum. Moreover, [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF] found some evidence suggesting that resonant wave interactions cannot explain the downscale cascade of energy of their turbulent simulations in systems close to geostrophy. In this study, however, we will not investigate this issue further.

In this study we calculate third-order structure functions from simulations of forced stratified rotating turbulence and compare the results from the simulations with the stratospheric data and the different analytical relations which we have reviewed in this introduction. In particular, we will focus on whether it is correct to make the approximation that are needed in order to derive (1.9) and (1.10). In addition to this, we will also investigate if our simulations are in agreement with another interesting observation made by Cho & Lindborg (2001) regarding the third-order horizontal velocity structure functions that are not parity invariant, that are the functions δu T δu T δu T and δu T δu L δu L , where δu T is the azimuthal velocity increment in the cyclonic direction. In the absence of system rotation parity invariance should hold and in this case these functions must be zero. In the presence of system rotation parity invariance is broken and there may be a cyclonic-anticyclonic asymmetry. Cho & Lindborg (2001) and Lindborg & Cho (2001) found that both these functions are positive in the lower stratosphere and display a clean r 2 -dependence in the interval 10 km < r < 1000 km. That these functions are positive means that cyclonic motions are dominant. However, no theoretical explanation can yet be offered for the r 2 -dependence.

The paper is organized as follows. In § 2, we briefly present the simulations. In § 3, we present the results for the third-order structure functions. Second-order statistics, such as energy spectra, have already been presented in [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. In this paper, we therefore focus on the third-order structure functions. In § 4, we end the paper with conclusions.

Simulations

The simulations are standard box simulations of the Boussinesq equations, using triply periodic boundary conditions and white noise (in time) forcing in a limited wavenumber band [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. The equations are formulated and solved for the three variables, q, a 1 and a 2 , where

q = - ∂u ∂y + ∂v ∂x + f N 2 ∂b ∂z , (2.1)
is the Charney potential vorticity [START_REF] Charney | Geostrophic turbulence[END_REF], and

a 1 = - f N ∂v ∂z + 1 N ∂b ∂x , (2.2a) a 2 = f N ∂u ∂z + 1 N ∂b ∂y (2.2b)
are two 'ageostrophic components' that measure the departure from thermal wind balance [START_REF] Vallis | Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation (Electronic Version)[END_REF]). The formulation using q, a 1 and a 2 is similar to the normal modes decomposition (see e.g. [START_REF] Bartello | Geostrophic adjustment and inverse cascades in rotating stratified turbulence[END_REF][START_REF] Embid | Low Froude number limiting dynamics for stably stratified flow with small or finite Rossby numbers[END_REF][START_REF] Smith | Generation of slow large scales in forced rotating stratified turbulence[END_REF]. Here u, v and w are the velocity components in the x, y and z directions, respectively, where x, y are the horizontal coordinates, z the vertical coordinate and b is the buoyancy. For the explicit form of the equations that are solved, the reader is referred to [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. The equations contain two dimensionless parameters, the Froude number and the Rossby number, which are defined as

Fr = U NL , Ro = U f L , (2.3a,b)
where U is a characteristic horizontal velocity, N the Brunt-Väisalä frequency, f the Coriolis parameter and L a characteristic horizontal length scale.

The system of equations is discretized in a triply periodic domain (L x , L y , L z ) = (2π, 2π, 2π), using QG scaling [START_REF] Charney | Geostrophic turbulence[END_REF], which means that the vertical coordinate has been stretched by a factor of N/f . The periodicity in all the three directions allows a Fourier representation along x, y and z; N x = N y = N z = 1024 modes were used in each direction. Pseudo-spectral methods are employed and the equations are advanced in time in spectral space. The nonlinear terms are evaluated in physical space and are discretized in time using a low-storage fourth-order Runge-Kutta scheme. In order to prevent aliasing errors, the 2/3-dealiasing rule was applied to the nonlinear terms. Linear terms are instead separately solved using an exact implicit procedure [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF].

Two different forcing schemes are adopted. In the first, which we call Q-forcing, the forcing is restricted to the q-equation. This is the same scheme as used in forcing or M-forcing, the forcing is applied to all three equations, with about twothirds of the energy injection into the ageostrophic modes, a 1 and a 2 . This scheme is similar to that used in a number of studies in stratified turbulence [START_REF] Brethouwer | Scaling analysis and simulation of strongly stratified turbulent flows[END_REF]Lindborg & Brethouwer 2007). In both cases energy is injected in a limited wavenumber band in the vertically stretched wavenumber space, k

∈ [k f -k/2, k f + k/2],
where k = 2. The characteristic length scale of the forcing is therefore L f = 2π/k f . As summarized in table 1 we have carried out three simulations at different Ro with large-scale (k f = 4) Q-forcing, three simulations at different Ro with large-scale M-forcing and one simulation with Q-forcing at an intermediate wavenumber (k f = 54). In all simulations the Froude number is equal to 0.01.

Energy is transferred both upscale and downscale, as compared to the forcing scale, with the ratio of energy cascading upscale and downscale mainly controlled by Ro [START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. It is worth noting that the ratio between the amount of energy cascading downscale and upscale is not a universal function and considerable differences are found between the Q-and the M-runs. Nevertheless, the same trend is observed, i.e. the ratio decreases with decreasing Ro. That M-runs have larger small-scale dissipation is also consistent with the analysis of [START_REF] Bartello | Geostrophic adjustment and inverse cascades in rotating stratified turbulence[END_REF] and the fact that ageostrophic energy (directly forced in the M-runs), by not being constrained by potential vorticity conservation, cascades directly to small scales. A large-scale linear drag as well as hyperviscosity (using a Laplacian to the power of four in a stretched grid) are introduced in the flow in order to allow for energy dissipation at large and small scales, respectively. After an initial transient, an approximate statistically stationary state develops from which statistics are collected.

Results

In figure 3 we plot the total spectral energy flux for the two runs Q2 and Q1. The curves are normalized with the total energy injection rate P. Part of the energy which is injected is going into an upscale energy cascade, which is seen as a negative flux at smaller wavenumbers than the forcing wavenumber. The remaining part is going into a downscale energy cascade which is seen as a positive flux at larger wavenumbers than the forcing wavenumber. As shown by [START_REF] Vallgren | A possible explanation of the atmospheric kinetic and potential energy spectra[END_REF], the fraction of energy which is going downscale increases with Ro which is also illustrated in figure 1. The total spectral flux can be divided into the flux of KE, Π K , and the flux of APE, Π P . There is a range where the fluxes are almost constant. We determine the spectral fluxes, Π K and Π P , by taking an average over k ∈ [15, 50]. In figure 4 we plot the function

- δu L δu • δu 2Π K r = - δu L δu L δu L + δu L δu T δu T 2Π K r (3.1)
for the six simulations with large-scale forcing. This function should be positive and equal to unity in a range where Π K is positive and where the relation (1.9) is satisfied. Indeed, all curves are positive in the relevant range of separation distances. If we used the curves in order to determine the sign of Π K we would therefore obtain the correct answer in all the cases. However, the curves generally overshoot unity. In all cases we would overestimate the magnitude of Π K by a factor of approximately 2 if we used relation (1.9) in a similar way as done by Cho & Lindborg (2001) for the third-order structure functions seen in figure 1.

In figure 5 we plot the function

- δu L δbδb 2N 2 Π P r (3.2)
for the same runs as in figure 4. For all curves there is a range where this function is reasonably close to unity, in agreement with the relation (1.10). For the simulation with Q-forcing at the lowest Rossby number, Ro = 0.05, this range is very narrow. However, at small separation distances the curve from this simulation collapses fairly well on the other curves. Thus, the results from the simulations support the interpretation of figure 2 as giving evidence of a positive or downscale energy flux of APE through the mesoscales, whose magnitude can be estimated by using the relation (1.10).

In figure 6 we plot ∇ • J K and ∇ • J P from run Q1. We also show the separate contributions from the horizontal divergence and the vertical divergence corresponding to the terms on the left-hand sides of (1.7) and (1.8). In addition, we show the energy exchange term 2 δwδb which appear with different signs on the right-hand sides of (1.5) and (1.6). In figure 6 normalized by Π P . The energy exchange term is small at small separations and in both figures there is a narrow range where the total normalized flux is equal to -4, confirming the theoretical relations (1.7) and (1.8). For the KE flux, the contribution from the vertical part of the divergence is half that from the horizontal divergence, but with opposite sign. This is the reason why the curves in figure 4 overshoot unity. On the other hand, for the APE flux, the contribution from the vertical part of the divergence is small and the total divergence can therefore be approximated by the horizontal part. This is the reason why the curves in figure 5 are all close to unity.

We find it somewhat unexpected that the vertical divergence makes such a significant contribution to the total divergence and has opposite sign to the horizontal divergence. Figure 7 shows the ratio (in absolute value) of the vertical and horizontal contributions to the total divergence for both Q-and M-runs. We have also carried out two additional simulations in order to check the convergence of the results with respect to the vertical resolution by decreasing the number of vertical points from 1024 to 512 and 256. In all cases, a plateau with values between 0.25 and 0.6 is found. The low-resolution runs show an increase of the ratio of only few percent.

On the other hand, the forcing scheme displays a stronger effect. The Q-runs shows a slightly larger dependence on Ro and diverge as the horizontal divergence vanishes. After this point, the ratio monotonically decreases and become small at large separation distance, indicating that the two-dimensional approximation would apply to a greater extent in this range of scales.

The third-order structure functions showed in figures 1 and 2 both display a switch of sign at a particular separation distance, r s . It is interesting to note that r s ≈ 400 km for the velocity structure function, in figure 1, as compared to r s ≈ 1000 km for the velocity-temperature-temperature structure function, displayed in figure 2. In figures 4 and 5 we can see that there is a corresponding switch of sign of the structure functions calculated from the simulations. However, for the third-order structure function including buoyancy in figure 5, r s is close to the forcing scale in all simulations, while for the velocity structure function in figure 4, r s moves towards smaller scales with decreasing Ro. This is particularly true for the simulations with Q-forcing.

In order to investigate this further we have decomposed the structure functions from the simulation into geostrophic and ageostrophic components, where the geostrophic component is the projection onto the subspace carrying only potential vorticity q and the ageostrophic component is the projection onto the complementary subspace [START_REF] Babin | On the asymptotic regimes and the strongly stratified limit of rotating Boussinesq equations[END_REF][START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]. In figure 8 we have plotted the velocity third-order structure function from runs Q1 (figure 8a) and M1 (figure 8b) including the separate contributions from the geostrophic and the ageostrophic components. In both cases, the contribution from the geostrophic component is positive and has an almost cubic dependence on r, while the contribution from the ageostrophic part is negative all the way up to the forcing scale and shows a linear dependence on r. This is particularly true for the Q1 run where the two components have comparable magnitude. In the M1 run on the other hand, the geostrophic part is much smaller and the ageostrophic part dominates throughout the whole range of scales. In figure 9 we see the corresponding curves for the buoyancy structure function. In run Q1, both the geostrophic contribution and the ageostrophic contribution are negative up to the forcing scale, where there is a switch of sign. The geostrophic contribution shows an approximate cubic dependence while the ageostrophic contribution shows a linear dependence on r. In the M1 run, on the other hand, the geostrophic part switches sign at intermediate scales, although its magnitude is several order of magnitude smaller than the ageostrophic part. In general, the Q-runs are characterized by a competition between geostrophic and ageostrophic dynamics whereas the M-runs are dominated by the dynamics of ageostrophic modes, the latter being due to wave dynamics and/or turbulent dynamics.

In the light of these results, we may interpret the switch of sign, at r s ≈ 400 km, of the velocity structure function in figure 1 as the transition scale between a regime which is dominated by geostrophic motions and a regime which is dominated by ageostrophic motions, while the corresponding switch of sign, at r s ≈ 1000 km, of the third-order structure function in figure 2, as a signature of the forcing scale in the stratosphere. As demonstrated in the recent analysis of GCM (general circulation model) data by [START_REF] Augier | A new formulation of the spectral energy budget of the atmosphere, with application to two high-resolution general circulation models[END_REF], the stratosphere is mainly forced by upwards travelling planetary waves carrying energy from the troposphere, where the main forcing mechanism is baroclinic instability.

Figure 4 shows that the scale at which δu L δu • δu switches sign varies significantly, depending on the forcing scheme and on Ro. Following the concept of a 'shadow' cascade [START_REF] Tung | The k -3 and k -5/3 energy spectrum of atmospheric turbulence: quasigeostrophic two-level model simulation[END_REF], [START_REF] Vallgren | A possible explanation of the atmospheric kinetic and potential energy spectra[END_REF] suggested that the transitional wavenumber k s can be identified as the scale at which the QG spectrum, E QG ∼ η 2/3 k -3 , with η being the potential enstrophy dissipation, and a stratified turbulence spectrum,

E ST ∼ ε 2/3 k -5/3 , attain comparable magnitude, k s ∼ √ η/ε, corresponding to a transitional scale r s = 2π/k s = 2π √ ε/η.
In figure 10(a) the transitional scale r s together with this prediction is shown for the Q-runs for which coexistence of a downscale enstrophy cascade and a downscale energy cascade is found (figure 10b). As figure 10(a) shows, the estimate of r s provides a fair agreement. In the Q-runs, the decrease of the small-scale dissipation with decreasing Ro leads to the observed reduction of r s .

In addition to the six simulation with forcing at k f = 4 we have also carried out a simulation with forcing at k f = 54. We have done this for two reasons. First, we would like to verify that our results are not critically dependent on the forcing wavenumber. Second, we would like to investigate to what extent the upscale cascade relation (3.3) (given below) for QG turbulence is valid. Quasi-geostrophic turbulence [START_REF] Charney | Geostrophic turbulence[END_REF]) is, in fact, very similar to two-dimensional turbulence. There are two quadratic invariants: total energy, which is the sum of KE and APE, and potential enstrophy, which is half the square of potential vorticity. Just as in two-dimensional turbulence, there are two possible cascades, a downscale cascade of potential enstrophy and an upscale cascade of energy. For the upscale energy cascade range of QG turbulence one can derive the relation [START_REF] Lindborg | Horizontal wavenumber spectra of vertical vorticity and horizontal divergence in the upper troposphere and lower stratosphere[END_REF])

δu L δu L δu L + δu L δu T δu T + δu L δbδb /N 2 = -2Π r, (3.3)
where Π is the total energy flux, which is negative since there is an upscale energy cascade. In the potential enstrophy cascade range there is no known exact third-order velocity structure function relation. However, for the enstrophy cascade range of twodimensional turbulence one can derive the exact relation [START_REF] Lindborg | Can the atmospheric kinetic energy spectrum be explained by two-dimensional turbulence?[END_REF])

δu L δu L δu L = δu L δu T δu T 1 8 Π ω r 3 , (3.4)
where Π ω is the enstrophy flux, which is positive. Cho & Lindborg (2001) fitted the measured third-order structure function in figure 1 to a positive cubic power law in a narrow range at separation distances of the order of 1000 km and interpreted this range as a possible enstrophy cascade range of two-dimensional turbulence. Using the relation (3.4) and the measured third-order structure function, Cho & Lindborg (2001) estimated the enstrophy flux as Π ω ≈ 2 × 10 -15 s -3 . This interpretation seems to be far more speculative than the interpretation of the linear negative range at mesoscales as evidence of downscale energy cascade. The range with a possible cubic dependence is clearly very narrow. At 1000 km scales, the assumptions that such an interpretation rest on also become questionable. One such assumption is the use of Taylor's hypothesis in the analysis of the aircraft data, that is that data taken at the two different points separated by distances of a 1000 km may be treated as simultaneous data. Another assumption is statistical homogeneity and isotropy in the horizontal plane. However, it is worth mentioning that the value which was estimated for Π ω is in close agreement with values estimated from GCMs [START_REF] Burgess | The troposphere-to-stratosphere transition in kinetic energy spectra and nonlinear spectral fluxes as seen in ECMWF analyses[END_REF].

In figure 11 we plot the spectral energy flux from two simulations with the same Ro but with the widely different forcing wavenumbers k f = 4 and k f = 54, respectively. In both simulations, roughly one tenth of the total energy injection goes into a downscale energy cascade while the remaining part goes into an upscale energy cascade. The reason why both the upscale and downscale energy fluxes are somewhat smaller in the k f = 4 simulation is that there is some direct action of the large-scale drag at the forcing scale in this simulation. In figure 12 we show the third-order structure functions for the k f = 54 simulation. The buoyancy structure function displays a change of sign at the forcing scale, with a negative range at smaller separations and a positive range at larger separation, while the velocity structure function displays a change of sign at a smaller scale than the forcing scale. This is consistent with the results of the simulations with k f = 4. As seen in figure 12, the QG relation (3.3) is quite well satisfied in the inverse cascade range. Cho & Lindborg (2001) and Lindborg & Cho (2001) calculated the transvere third-order velocity structure functions δu T δu T δu T and δu T δu L δu L . In these studies the transverse component was taken as positive in the anticyclonic direction, that is the clockwise azimuthal component on the northern hemisphere and the anti-clockwise component on the southern hemisphere. In the following we will do it the other way around and take the transverse component as positive when it is in the cyclonic direction. With this convection, Cho & Lindborg (2001) found that both these functions are positive and show a quadratic dependence on separation distance in the lower stratosphere. Thus, there is a preference for cyclonic motions in the lower stratosphere.

In figure 13 we have plotted the sum of the two functions, δu T δu T δu T + δu T δu L δu L , normalized by E 3/2 K , where E K is the mean KE. In all simulations this function is positive and shows an approximate quadratic dependence on r, in agreement with the results from the stratosphere. There is no clear trend with regard to Ro and there is no significant difference between the simulations with Q-forcing and M-forcing. The magnitude of the curves from the simulations with Q-forcing is smaller, but this is due to the normalization, since E K is larger in these simulations. In figure 14 we have plotted the skewness δu 3 T / δu 2 T 3/2 from the simulations and also inserted a curve which is calculated based on the results given by Lindborg & Cho (2001). For the atmospheric data, we have normalized r using L f = 1000 km. As seen in figure 14, the skewness corresponding to the stratospheric measurements is a little bit larger than unity, indicating that there is a strong dominance of cyclonic motions, while it is lower, ∼0.1-0.5, in the simulations. We can offer no theoretical explanation of these results, but only point out that there is both an interesting similarity between the observations and the simulations, in that there is a preference for cyclonic motions and that the curves show a quadratic dependence on r, and a striking dissimilarity, in that the magnitude of the skewness is considerably larger for the stratospheric data compared to the simulations.

Conclusions

Cho & Lindborg (2001) measured the third-order velocity structure functions using aircraft data in the lower stratosphere. They found that δv L δv L δv L + δv L δv T δv T has a negative linear dependence on r in the mesoscale range and interpreted this as a sign of a downscale cascade of KE. They used the two-dimensional relation 4; black line: measurements in the atmosphere (evaluated from Cho & Lindborg 2001). The separation distance for the atmospheric data has been normalized using L f = 1000 km.

(1.9) to estimate the spectral KE flux, Π , through the scales. In this paper we have suggested that the negative linear third-order velocity-temperature-temperature structure function measured by [START_REF] Lindborg | Determining the cascade of passive scalar variance in the lower stratosphere[END_REF] may be interpreted as a sign of a downscale cascade of APE. Using the two-dimensional flux relation (1.10), the spectral energy flux of APE, Π P , may be estimated. In order to investigate to what extent the interpretation of the measurements can be valid we have carried out idealized simulations of rotating and stratified turbulence at different Rossby numbers, using two different forcing schemes. In all simulations, the method used by Cho & Lindborg (2001) gives a correct prediction of the sign of Π K . However, the magnitude is overestimated by a factor of two. The reason for this is that the term including the vertical derivative in (1.7) cannot be neglected. For the interpretation of the third-order structure function including buoyancy, the use of the two-dimensional flux relation (1.10) provides a much more accurate estimate of the downscale energy flux.

We have also calculated the third-order structure functions which are not parity invariant and found that these scale as r 2 , similar to what was found in the stratosphere, with a strong preference for cyclonic motions. However, the skewness of the cyclonic velocity increment, δu T , is lower in the simulations than in the stratosphere. There are a number of studies of homogeneous rotating turbulence showing a preference for cyclonic motions (e.g. [START_REF] Hopfinger | Vortices in rotating fluids[END_REF][START_REF] Bartello | Coherent structures in rotating three-dimensional turbulence[END_REF]. In these studies, the skewness of vertical vorticity, ω z , is generally used as a measure of the degree of symmetry breaking. In a field with an energy spectrum which is more shallow that k -3 , vorticity is a small-scale quantity. Therefore, the skewness of ω z may be influenced by the small scales and may thus be Reynolds-number dependent [START_REF] Moisy | Decay laws, anisotropy and cycloneanticyclone asymmetry in decaying rotating turbulence[END_REF]. On the other hand, the skewness of δu T is a measure of the symmetry breaking as a function of scale and is therefore more informative. We suggest that this measure should be used in order to study cyclonic/anticyclonic symmetry breaking.

The most important conclusion of our study is that it seems justified to interpret the measured third-order structure functions in the lower stratosphere as giving evidence of a downscale energy cascade of KE and APE, which is similar to what is seen in the simulations. In particular, we think that the very clean observational result for the velocity-temperature-temperature structure function, in figure 2, together with the corresponding result from the simulations, in figure 5, provides support for interpreting this as a result of the downscale APE flux, allowing us to find the estimate of Π P ≈ 2 × 10 -5 m 2 s -3 . Had the same third-order structure functions as measured by Cho & Lindborg (2001) displayed a positive sign or a dependence on r far from linear in our simulations, then there would be no justification for interpreting the measurements of Cho & Lindborg (2001) as evidence of a downscale energy cascade. Cho & Lindborg (2001) also measured the third-order structure functions in the upper troposphere. However, in this case it was difficult to obtain any consistent and well-converged results to be interpreted in the light of turbulence theory. There are probably two reasons why the results from the upper troposphere were less clean. First, the Brunt-Väisälä frequency shows a rather significant variation with height in the upper troposphere [START_REF] Charney | Propagation of planetary-scale disturbances from the lower into the upper atmosphere[END_REF], which may invalidate the assumption of statistical homogeneity in the vertical direction. Second, the release of latent heat plays a more important role in the upper tropospheric dynamics than in the stratospheric dynamics. In addition to this, one may also speculate that the difference of the forcing mechanism leads to some difference in the cascade dynamics. The troposphere is mainly forced by baroclinic instability, while the stratosphere is forced by upwards propagating planetary waves. The lower stratosphere seems to be an almost ideal test bed for theories on rotating and stratified turbulence. Generally, the Brunt-Väisälä frequency shows a small variation with height, so that statistical homogeneity is fairly well satisfied in the vertical direction and the release of latent heat is of minor importance. It is our hope that the present study will inspire more observational work on rotating and stratified turbulence in the lower stratosphere.

FIGURE 1 .

 1 FIGURE 1. Sum of the diagonal third-order structure functions. Straight lines indicate the power laws -r and r 3 . Circles and crosses indicate negative and positive values, respectively. (Reproduced from Cho & Lindborg 2001.)

FIGURE 2 .

 2 FIGURE 2. Plot of δu T δTδT /(2r) versus r. Solid line: best linear fit, 20 km < r < 400 km. Circles and crosses indicate negative and positive values, respectively. (Reproduced from Lindborg & Cho 2000.)

  FIGURE 3. (Colour online) Total energy flux normalized with the total energy input P for --, Q1 (Ro = 0.1) and -• -, Q2 (Ro = 0.2). The forcing scale is at k f = 4.

FIGURE 4 .FIGURE 5 .

 45 FIGURE 4. Plot ofδu L δu • δu /2Π K r for the different runs. Solid lines: Q-forcing; dashed lines: M-forcing. Green Ro = 0.05; red Ro = 0.1; blue Ro = 0.2.

  FIGURE 6. (Colour online) Decomposition of the divergence of (a) J K = δvδu • δu and (b) J P = δvδbδb /N 2 into the horizontal contribution (----) and the vertical contribution (-• -) for run Q1. Also shown are the sum of the horizontal and vertical contributions (--) and the conversion from APE to KE 2 (• • • • • • • •). Grey lines indicate the theoretical prediction -4.

FIGURE 7 .FIGURE 8 .

 78 FIGURE 7. Ratio between the vertical, ∂/∂z δwδu • δu , and the horizontal contributions (1/r h )(∂/∂r h ) r h δuδu • δu to the total divergence as a function of separation distance r. Coloured lines as in figure 4. Black lines refer to simulations in which the vertical resolution has been decreased to 512 (-• -) and 256 (• • • • • • • •).

FIGURE 9 .

 9 FIGURE 9. Decomposition of δu L δbδb (blue) into its geostrophic component (red) and ageostrophic component (green). (a) Q1 and (b) M1. Solid lines: positive values; Dashed lines: negative values. Grey lines as in figure 8: solid, 2Π P r; dashed, ∼ r 3 .

  FIGURE 10. (Colour online) (a) Transition scale r s at which δu L δu • δu changes sign. , Q-runs; ×, prediction of the transitional scale r s = 2π √ ε/η (Vallgren et al. 2011). (b) Flux of energy (--) and enstrophy (----) for the Q05 run. Vertical dotted line represents the transitional wavenumber k s = 2π/r s ≈ 50, i.e. k s ≈ 10k f .

FIGURE 11

 11 FIGURE 11. (Colour online) Total energy flux normalized with the total energy input P for Ro = 0.1. --, Run Q1 (k f = 4); -• -, run Q1-SS (k f = 54).

FIGURE 12 .FIGURE 13 .

 1213 FIGURE 12. Plot of J K (red), J P (green) and J K + J P (blue) for run Q1-SS (k f = 54): --, positive values; ----, negative values. Solid grey line, (3.3); dashed grey line to aid the reader, ∼ r 3 .

TABLE 1 .

 1 Summary of the simulations. The last column shows the ratio between small-scale dissipation, ε, and energy injection, P, in the quasi-stationary state.

	Run	Ro	k f	ε/P
	Q2 Q1 Q05 Q1-SS 0.1 54 0.2 4 0.1 4 0.05 4 5.9 × 10 -3 0.27 0.08 0.08
	M2 M1 M05	0.2 0.1 0.05 4 4 4	0.88 0.78 0.69

[START_REF] Vallgren | A possible explanation of the atmospheric kinetic and potential energy spectra[END_REF]

and

[START_REF] Deusebio | The route to dissipation in strongly stratified and rotating flows[END_REF]

. In the second scheme, which we call mixed

Acknowledgement

Computer time was provided by Swedish National Infrastructure for Computing (SNIC) which is gratefully acknowledged.


Sarmad Munir
email: sarmad.munir@itk.ntnu.no

Morten Hovd
email: morten.hovd@itk.ntnu.no

Sorin Olaru
email: sorin.olaru@centralesupelec.fr

Low complexity constrained control using higher degree Lyapunov functions

Keywords: Contractive Sets, Constrained Control, explicit MPC, Ellipsoidal Contractive sets

Explicit Model Predictive Control often has a complex solution in terms of the number of regions required to define the solution and the corresponding memory requirement to represent the solution in the online implementation. An alternative approach to constrained control is based on the use of controlled contractive sets. However, polytopic controlled contractive sets may themselves be relatively complex, leading to a complex explicit solution, and the polytopic structure can limit the size of the controlled contractive set. This paper develops a method to obtain a larger controlled contractive set by allowing higher order functions in the definition of the contractive set, and explores the use of such higer-order contractive sets in controller design leading to a low complexity explicit control formulation.

Introduction

The ability to capture operational constraints is of vital importance in controller design for real-life applications. It is reasonable to state that the ability to handle constraints in a transparent way is what sets the industrially very successful Model Predictive Control (MPC) [START_REF] Qin | A Survey of Industrial Model Predictive Control Technology[END_REF] apart from the theoretically elegant -but less industrially successful -LQG control. Standard MPC solves an optimization problem online, but due to the computational complexity of MPC it is limited to the systems which are not safety critical (due to the use of complex and thus error prone optimization software), have sufficiently slow dynamics, and/or can afford high performance computational hardware [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF]. Explicit MPC [START_REF] Bemporad | The explicit linear quadratic regulator for constrained systems[END_REF] to some degree resolves this problem and allows the use of low-complexity computing code in the online implementation. Unfortunately, the explicit solution to standard MPC problems often has a highly complex solution, and even in cases when the explicit solution can be found in acceptable time the implementation of the solution on the online control hardware may require excessive memory. Low complexity constrained control with modest computational complexity, small memory requirements and simple, thus verifiable code in the online implementation is therefore desired.

One approach to such low complexity constrained control is based on the use of a controlled contractive set. The complexity of the solution will then depend on the complexity of the contractive set. Therefore, obtaining a controlled contractive set of low complexity is essential for this approach to formulate low complexity explicit constrained control. A maximal polyhedral controlled contractive set with a given contraction factor can be obtained by the iterative procedure described in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF]. However, the complexity of the contractive set thus obtained may be very high. A non-iterative procedure for obtaining a contractive set of low complexity is proposed in [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF]. The approach is not applicable to systems with identical modes in series (corresponding to a non-diagonalizable A-matrix in the system's state space representation). Furthermore, the contractive set obtained in [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF] is of fixed complexity, which does not allow trading off the complexity against the size of the contractive set. An optimization based technique has been proposed in [START_REF] Munir | Controlled contractive sets for low-complexity constrained control[END_REF] which allows the trading off complexity versus the size of the set. A solution to the optimization problem in [START_REF] Munir | Controlled contractive sets for low-complexity constrained control[END_REF] not only reduces the on-line computational complexity of the resulting constrained control, but also ensures significant reduction in the memory required to store the explicit solutions. However the method explained in [START_REF] Munir | Controlled contractive sets for low-complexity constrained control[END_REF] is highly non-convex, which makes it difficult to use for finding sufficiently large contractive sets for higher dimensional systems. Alternatively, ellipsoidal contractive sets with corresponding linear control laws can be computed, but the measure of these sets is limited by the linear structure of the control law and the inherent conservatism of the corresponding quadratic Lyapunov function. This paper proposes a method to obtain an enlarged contractive set by defining the contractive set using a function of variable degree (a degree which is assumed to be greater or equal to 2, thus including the quadratic forms as a particular case), and also allowing for higher order control laws. Note that the function defining the controlled contractive set can be interpreted as controlled Lyapunov function for the closed loop system.

In Section 2, the controller design using controlled contractive sets is presented, along with the formulation for finding the largest ellipsoidal controlled contractive set fulfilling state and input constraints. Section 3 describes a controller design which leads to the determination of larger contractive sets. The controller design is inspired by the results in [START_REF] Oishi | Direct Design of a Polynomial Model Predictive Controller[END_REF], but unlike the respective work, the controller will be defined using only two regions. The method described in Section 3 is applied to illustrative examples and the results are described in Section 4, which is followed by a discussion and conclusions in Section 5.

Contractive Sets

Consider the constrained control of the linear discrete time system:

x k+1 = Ax k + Bu k (1)
with x k ∈ nx , u k ∈ nu representing the current state and input, respectively, while x k+1 is the next time step state. The system is subject to input constraints

U = {u k |H u u k ≤ 1}, with H u ∈ np u ×nu .
Definition 1 Given a function V : nx → , the level set of V (x) for a scalar α is the set S α = {x|V (x) ≤ α}.

Proposition 1 Consider a function V (x) : nx → satisfying the following properties:

A1 positive definite, with V (0) = 0, A2 continuous, A3 radially unbounded, i.e., V (x) → ∞ as x → ∞.

Then

(1) All level sets S α exist and are bounded for all 0 ≤ α < ∞. (2) If β < α, S β ⊂ S α .

Proof:

From A1 it follows that the level sets S α = ∅ if α < 0. Claim (1) follows directly from A1, A2 and A3. For claim [START_REF] Blanchini | Set Theoretic Methods in Control -Second Edition[END_REF] we note that S β ⊆ S α is a consequence of Definition 1. Next, consider two points x 1 and x 2 with V (x 1 ) = V (x 2 ) + δ for some δ > 0. Then by applying the Mean Value Theorem, continuity of V (x) implies that the points x 1 and x 2 must be separated by some nonzero distance. Hence, we get strict inclusion,

S β ⊂ S α if β < α. 2
Definition 2 Consider a continuous and radially unbounded function V : nx → ≥0 . A level set S α is controlled γ-contractive with respect to (1) for a given γ ∈ (0, 1), if ∀x k ∈ S α , ∃u k ∈ U such that x k+1 ∈ S γα .

The functions V (x) fulfilling the assumptions of Proposition 1 are natural ingredients in control designs enforcing contractiveness properties, as for example in the low complexity optimization based formulation min

u k ,x k+1 1 2 x T k+1 Qx k+1 + 1 2 u T k Ru k (2a) subject to x k+1 = Ax k + Bu k (2b) H u u k ≤ 1 (2c) V (x k+1 ) ≤ γV (x k ) (2d)
where Q and R represent the state and input weights.

Consider next the bounded state constraints

x k ∈ X with X = {x k |H x x k ≤ 1} where H x ∈ px×nx .
Proposition 2 Let V (x) be a function fulfilling assumptions A1 -A3 of Proposition 1, and let V (x) = α, Then, if

(1) the corresponding level set S α is controlled γcontractive, and (2) S α ⊆ X the control action obtained as a solution of (2) guarantees an exponentially stability of the closed loop which in addition fulfills input and state constraints over S α .

Proof: Follows directly from Proposition 1 and Definition 2. 2

As a result of Proposition 2, the function V (x) is a Lyapunov function for the system (1) inside the set S ᾱ, where ᾱ = max α such that S α ⊆ X .

In [START_REF] Hovd | Low Complexity Constraint Control Using Contractive Sets[END_REF] and [START_REF] Munir | Controlled contractive sets for low-complexity constrained control[END_REF], a controller based on (2) with polytopic controlled contractive sets S = {x k |F x k ≤ 1} were studied based on a piecewise linear function

V (x k ) = max{F x k } (3)
Using the function specified as in (3), the optimization (2) becomes a standard quadratic program, which may be solved parametrically with x k and V (x) as parameters. This is done by imposing a virtual parameter α k = V (x k ) before solving the optimization in (2) at time k. The constraint in (2d) then simply becomes

F (Ax k + Bu k ) ≤ γα k .
As the total number of constraints and the number of degrees of freedom are typically quite modest in (2) compared to a classical MPC problem utilizing a longer prediction horizon, the parametric solution is also of modest complexity. However, this approach suffers from the drawbacks described in the Introduction, and this paper therefore focuses on allowing more general types of function V (x), to obtain a larger operating region with modest online computational complexity and memory requirement for the control.

In the developments below, two ellipsoidal controlled contractive sets will be important as terms of comparison:

• The set Ω = {x ∈ nx |x T P -1 x ≤ 1}, the largest controlled γ-contractive set that can be obtained using linear state feedback. • The set Ω uc = {x ∈ nx |x T P -1 uc x ≤ 1}, the ellipsoidal set where γ-contractiveness is achieved with the linear state feedback u k = K uc x k .

Constraints in both states and inputs are accounted for in the calculation of both Ω and Ω uc . These sets can be calculated using well known techniques based on Linear Matrix Inequalities, see, e.g., [START_REF] Blanchini | Set Theoretic Methods in Control -Second Edition[END_REF] or [START_REF] Ngyuyen | Constrained control of uncertain systems: An interpolation baszed approach[END_REF] for details.

While the set Ω uc can be found for any given controller 1 K uc , for the subsequent use in this paper it will be considered to be the unconstrained solution to [START_REF] Blanchini | Set Theoretic Methods in Control -Second Edition[END_REF], see the inverse optimality arguments in [START_REF] Nguyen | An inverse optimality argument to improve robustness in constrained control[END_REF] for the choice of weights Q and R. When ignoring the input and contractivity constraints, (2) yields the controller

u k = -(R + B T QB) -1 B T QA Kuc x k
For notational convenience in the following, we will define P 1 = P -1 and P 0 = P -1 uc .

3 Controller Design for Higher Order Contractive Set

Approximate optimization problem solution

Despite the computational advantages of the ellipsoidal sets recalled above, the associated linear feedback and the quadratic structure of the Lyapunov function limit the volume of the contractive set. Larger contractive sets can be obtained by increasing the complexity of the Lyapunov function. A natural approach is to focus on polynomial forms of higher order. By relaxing the control structure to allow the input to have higher order dependency on the state, further relaxations on the contractive set can be obtained.

Using (2b) to eliminate x k+1 from the optimization formulation, (2) may be reformulated as Problem P u :

min u k 1 2 u T k Hu k + x T k F u k (4a) subject to H u u k ≤ 1 (4b) V (x k+1 ) ≤ γV (x k ) (4c) where H = (B T QB + R), F = A T QB.

The Lagrangian function for Problem

P u is L(u k ) = 1 2 u T k Hu k + x T k F u k +λ T u (H u u k -1) + λ q (V (x k+1 ) -γV (x k )) (5)
The corresponding KKT conditions are

Hu k + F T x k + H T u λ u + V 1 k+1 (u k )λ q = 0 (6a) H u u k -1 ≤ 0 (6b) V (x k+1 ) -γV (x k ) ≤ 0 (6c) λ u ≥ 0 (6d) λ q ≥ 0 (6e) λ T u (H u u k -1) + λ q (V (x k+1 ) -γV (x k )) = 0 (6f) where V 1 k+1 (u k ) denotes the first derivative of V (x k+1 ) with respect to u k . Next, consider the optimizaton problem P c : min c,u k ,λu,λq c (7a)
subject to (6a) -(6e) and

-λ T u (H u u k -1) -λ q (V (x k+1 ) -γV (x k )) ≤ c (7b)
Denote the solution to (7) by c * , u * k , λ * u , λ * q . Observe that if c * = 0, u * k is also the optimal solution to P u . Clearly, c * < 0 is not possible. The constraints of Problem P c includes the constraints of Problem P u , and hence a feasible solution to P c is also a feasible solution to P u with u * k (obviously fulfilling the constraints of problem P c ) being a suboptimal solution to P u .

Define by

J(u * k (c)) = 1 2 (u * k ) T Hu * k + x T k F u *
k the value of the objective function of Problem P u evaluated with the input u * k from the solution to Problem P c . Correspondingly, let J(u * k (0)) denote the (optimal) value function for Problem P u .

Lemma 3 Consider an optimal solution (c * , u * k , λ * u , λ * q ) to P c , with c * > 0. Then u * k is a suboptimal solution to P u , with J(u * k (c)) -J(u * k (0)) < c.

Proof:

This proof follows the approach in [START_REF] Oishi | Direct Design of a Polynomial Model Predictive Controller[END_REF]. For any feasible

u k 1 2 u T k Hu k + x T k F u k ≥ 1 2 u T k Hu k + x T k F u k + λ * u λ * q T H u u k -1 V (x k+1 ) -γV (x k ) = M (u k ) (8)
Next, minimize both sides subject to constraints (6b) and (6c). Thus

J(u * k (0)) ≥ min (6b),(6c) M (u k ) ≥ min u k ∈ nu M (u k ) (9)
The function M (u k ) can be recognized as the Lagrangian function or P u , with fixed λ u = λ * u and λ q = λ * q . The unconstrained minimization of M (u k ) yields (6a), again with λ u = λ * u and λ q = λ * q . One therefore finds that the optimal value for the unconstrained minimization of

M (u k ) yields u k = u * k . Thus, J(u * k (0)) ≥ M (u * k )
Multiply the inequality above with -1 and add J(u * k (c)) to both sides to obtain

J(u * k (c)) -J(u * k (0)) ≤ c * (10) 2

Problem reformulation

In the following, the approximate solution to the optimization problem P u is sought as a function of the present state vector x k . However, for u k expressed as a polynomial function of x k , V (x k+1 ) will have a higher order dependence on the polynomial coefficients of u k (x k ), which will cause problems in the resulting controller design formulation. To circumvent this problem, attention is refocused on the equivalent formulation (2), the KKT conditions of which are

Ru k -B T λ e + H T u λ u = 0 (11a) Qx k+1 + λ e + λ q ∇V (x k+1 ) = 0 (11b) x k+1 -Ax k -Bu k = 0 (11c) H u u k -1 ≤ 0 (11d) V (x k+1 ) -γV (x k ) ≤ 0 (11e) λ u ≥ 0 (11f) λ q ≥ 0 (11g) -λ T u (H u u k -1) -λ q (V (x k+1 ) -γV (x k )) = 0 (11h)
where the operator ∇ ∆ = d dx k+1 (.). As before, the complementarity constraints are relaxed, and solutions with a commensurate relaxation are sought, in which the input is expressed as a function of the present state x k . However, as a novelty of the present approach, instead of considering the model equations as a separate constraint, this will be added as an extra term in all other constraints involving both x k and x k+1 . This yields the formulation min c (12a) subject to constraints (11d), (11f), (11g), and

Ru k -B T λ e +H T u λ u +µ T 1 (x k+1 -Ax k -Bu k ) = 0 (12b) Qx k+1 + λ e + λ q ∇V (x k+1 ) + µ T 2 (x k+1 -Ax k -Bu k ) = 0 (12c) V (x k+1 ) -γV (x k ) + µ T 3 (x k+1 -Ax k -Bu k ) ≤ 0 (12d) -λ T u (H u u k -1) -λ q (V (x k+1 ) -γV (x k )) + µ T 4 (x k+1 -Ax k -Bu k ) ≤ c (12e)
The multipliers µ 1 , µ 2 , µ 3 and µ 4 as well as λ e are polynomial functions of x k and x k+1 with no positivity constraint. The multipliers λ q and λ u are also polynomials in x k and x k+1 , but have the positivity constraints as explicitly stated above. We retain from the control framework that it is important to evaluate the input based on information available online, and hence u k can only be designed as a function of x k . The degrees of freedom in this optimization are the polynomial coefficients of the multipliers, the coefficients of the input considered as a polynomial feedback function, and finally the coefficients of the polynomial Lyapunov function V (x k ). This reformulation can be interpreted as lifting the original problem formulation to a higher dimensional space with both x k and x k+1 as independent variables, while the terms involving the µ i ensure that the equations in which they are inserted hold on the manifold of system trajectories (where the term inside the parentheses is identically zero). The system equations (1) are fulfilled by the physics of the system -and thus need not be enforced by the control. In contrast, the controller design has to ensure that equalities (11a) and (11b) are fulfilled along the trajectories of the system, in order to ensure an approximately optimal control.

The controller design attempts to find a feasible approximately optimal solution to (2) for the set

S = {x|V (x) ≤ 1}
The design procedure is initialized with S = Ω and V (x) = x T P 1 x, and thereafter the set S is iteratively enlarged while allowing for higher order V (x). The optimal controller inside Ω uc is already known to be K uc . Note that Ω uc ⊆ S. The controller design can therefore be divided in two parts, one is the unconstrained controller for Ω uc and the other controller is designed for a region inside S but outside the ellipsoid Ω uc . This region can be defined by the set

S C = {x k |p(x k ) > 0} (13a)
where

p(x k ) = -(1 -V (x k ))(1 -x T k P 0 x k ) (13b)
Using the S-procedure, the constraints defined above can be enforced in the region where p(x k ) is positive. In addition, it is natural to allow for a larger relaxation of the complementarity constraints (and thus larger absolute distance to optimum) when the state is far from the origin. Therefore, the complementarity constraints are relaxed by a factor cx T k x k instead of relaxing only by c. This yields the optimization formulation min c (14a)

Ru k -B T λ e +H T u λ u +µ T 1 (x k+1 -Ax k -Bu k ) = 0 (14b) Qx k+1 + λ e + λ q ∇V (x k+1 ) + µ T 2 (x k+1 -Ax k -Bu k ) = 0 (14c) H u u k -1 + s 1 (x k )p(x k ) ≤ 0 (14d) V (x k+1 ) -γV (x k )+µ T 3 (x k+1 -Ax k -Bu k ) + s 2 (x k , x k+1 )p(x k ) ≤ 0 (14e) λ u -s 3 (x k )p(x k ) ≥ 0 (14f) λ q -s 4 (x k )p(x k ) ≥ 0 (14g) -λ T u (H u u k -1) -λ q (V (x k+1 ) -γV (x k ))+ µ T 4 (x k+1 -Ax k -Bu k ) + s 5 (x k , x k+1 )p(x k ) ≤ cx T k x k (14h)
Here s 1 (x k ), s 2 (x k , x k+1 ),s 3 (x k ), s 4 (x k ) and s 5 (x k , x k+1 ) are SOS polynomials.

It was explained above that the set S is increased iteratively. Note that S j-1 ⊆ S j provided:

(1 -V j (x k )) -s 6 (x k )(1 -V j-1 (x k )) ≥ 0 (15)
It is not yet guaranteed that the set S fulfills the state constraints. The state constraints are given by H x x k ≤ 1 and thus may be regarded as the intersection of the n c unbounded sets C r = {x k |H x,r x k ≤ 1}, where H x,r is the r th row of H x . Hence the state constraints become:

(1-H x,r x k )-σ r (x k )(1-V j (x k )) ≥ 0, r = 1, . . . p x ( 16
)
where s 6 (x k ) and σ r (x k ) are SOS polynomials. Therefore, an approximate solution to the optimization problem described in (2) can be found by solving

min c (17a) subject to Ru k -B T λ e +H T u λ u +µ T 1 (x k+1 -Ax k -Bu k ) = 0 (17b) Qx k+1 + λ e + λ q ∇V (x k+1 ) + µ T 2 (x k+1 -Ax k -Bu k ) = 0 (17c) H u u k -1 + s 1 (x k )p(x k ) ≤ 0 (17d) V (x k+1 ) -γV (x k )+µ T 3 (x k+1 -Ax k -Bu k ) + s 2 (x k , x k+1 )p(x k ) ≤ 0 (17e) λ u -s 3 (x k )p(x k ) ≥ 0 (17f) λ q -s 4 (x k )p(x k ) ≥ 0 (17g) -λ T u (H u u k -1) -λ q (V (x k+1 ) -γV (x k ))+ µ T 4 (x k+1 -Ax k -Bu k ) + s 5 (x k , x k+1 )p(x k ) ≤ cx T k x k (17h) (1 -V j (x k )) -s 6 (x k )(1 -V j-1 (x k )) ≥ 0 (17i) (1 -H x,r x k ) -σ r (x k )(1 -V j (x k )) ≥ 0, r = 1, . . . n c (17j)

Solving the Sums-of-Squares problem

The optimization problem described in (17) can be solved by using sum-of-squares (SOS) programming.

The SOS method has received a lot of attention since the PhD thesis of Parrilo [START_REF] Parrilo | Structured Semidefinite Programs and Semialgebraic Geometry Methods in Robustness and Optimization[END_REF]. The SOS technique generalizes the known algorithmic tool in Linear Matrix Inequalities (LMIs) for which there exist many efficient solvers. The SOS method employs the similar techniques as LMI problems, but all problems are formulated in terms of the polynomials or polynomial matrices [START_REF] Papachristodoulou | A Tutorial on Sum of Squares Techniques for Systems Analysis[END_REF].

In the problem described in (17), all constraints are of SOS type except the constraints (17b) and (17c), which gives linear constraints on the parameters of the polynomials u k (x k ), λ e (x k , x k+1 ), λ q (x k , x k+1 ) and λ u (x k ). The polynomials s 1 ,s 2 ,s 3 ,s 4 ,s 5 , s 6 and σ r are all SOS polynomials in x k , as is also V (x k ). It can be noticed that the problem is bilinear in the coefficients of s i , σ r and V (x k ). It is also bilinear with respect to coefficients of λ u and u k . Instead of attempting to solve the bilinear problem formulation directly, this paper instead uses the common approach of iteratively solving linear subproblems:

• First λ u , λ q , λ e , s i and σ r are optimized with given V (x k ) and u k . • Then V (x k ) and u k are optimized, with λ u , λ q , s i and σ r from the step above.

The problem (17) makes sure that the set S j-1 ⊆ S j but it does not guarantee that S j-1 ⊂ S j and it may result in S j-1 = S j for any number of iterations. To enforce that the set S j is larger than the set S j-1 , we select the points outside S j-1 and ensure that they are included in S j , which results in a larger set S j . After a certain number of iterations we may not be able to find a set S j greater than S j-1 . Increasing the degree of V (x k ) and/or u(x k ), one may then be attempted in order to increase the size of S j further.

Algorithm

This section provides a comprehensive procedure to obtain a controlled contractive set. We start with the linear system and obtain the largest ellipsoidal contractive set Ω and the largest unconstrained ellipsoidal contractive Ω uc set for that system as discussed in section 2. These sets always exist in the linear system framework. We want to find a contractive set S which is larger than Ω. Of interest is only the part of this contractive set where the function p(x k ) is positive. This controller will ensure that all the system trajectories contract towards the origin at least by factor γ. Once the trajectories are inside Ω uc , the controller K uc will take over (representing in fact a switched control with a state-based switching rule).

We can choose the degree x deg of V (x k ) (which also defines the set S) and the degree u deg of the control law u k as per system requirements. Select num number of uniformly distributed points on the boundary of the ellipsoid Ω and check which point undergoes most contraction. In this way we select a point where the contraction constraint is farthest from being active. The idea is to push this point outwards and ensure that S contains the new scaled point, in this way S will be larger than Ω. As discussed in the section 3.3, obtaining S is an iterative procedure. In the first iteration, the problem formulated in (17) is solved by keeping S = Ω and u k as the corresponding control law for Ω. By solving the problem, the optimized values for λ u , λ q , λ e , µ, s i , and σ r can be obtained. These are now kept fixed and(17) can be updated by solving for V (x k ) and u k . An additional constraint is added to (17), which if a valid solution is found ensures that the point outside Ω are included in the set S. Now we can take the new set S as our starting set and repeat the whole procedure to obtain even larger set. This process should continue until the contractive set cannot be enlarged any further. In that case the iteration should stop, or either x deg or u deg (or both) should be increased.

Set j = j + 1 4:
Find boundary points of the set S j-1 by solving for V (x k ) == 1 along rays in directions defined by vectors from the origin to points uniformly distributed on the n x -dimensional unit sphere.

5:

Check which point contracts most by applying control law u k . Select that point point.

6:

Solve the SOS problem (17) by keeping V (x k ) and u k fixed. Optimize for λ u , λ q , λ e , s i and σ r .

7:

Set a small positive number increment, which will be used to specify a point outside S j-1 .

8:

Find a point outside S j-1 by adding increment × point to point i.e. point new = point+increment× point.

9:

Keep λ u , λ q , s i and σ r fixed, optimize (17) using V (x k ) and u k , while ensuring that the point point new is included in S j . Update V (x k ), u k , λ u , λ q and s i .

12: end if 13: end while Let the controller obtained from algorithm 1 be denoted by u k = c u (x k ), then the explicit controller for the overall problem becomes:

u k = K uc x k , if x k ∈ Ω uc c u (x k ), if x k ∈ S C (18)

Examples

The application of the method described in section 3 is illustrated in this section. The SOS problems are formulated and solved using YALMIP [START_REF] Löfberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF].

Example 1

Consider a linear system whose state representation is given as:

x k+1 = 1 0.2 0 1 x k + 0.22 0.2 u k (19)
The input constraints are given as -2 ≤ u k ≤ 2 and each element j of the state vector should fulfill the constraints The state trajectories for the dark gray set are shown in figure 2. The method described in section 3 is used to design the controller for the region (dark gray), which is outside the ellipsoidal (light gray) set, whereas inside the ellipsoidal set the unconstrained controller K u c is used. The ellipsoidal set is the set corresponding to the unconstrained controller, and the unconstrained (linear) control law is applied once state is inside the ellipsoidal set. The value of c, the measuring distance from the optimum comes out to be 4.54 • 10 -12 . In Figure 3, contractive sets with variable input degree are shown. All sets use a Lyapunov function of degree 8. The white set is found with a controller of degree 3, while the dark and light gray sets correspond to the controllers of degree 2 and 1 respectively.

Example 2

Consider the following system with oscillatory modes (corresponding to a complex conjugate pair of eigenvalues).

x

k+1 = 0.3 -0.2 0.1 0.1 x k + 0.2 0.4 0.4 0.5 u k (20)
The input and the state constraints are given as -2 ≤ u k,i ≤ 2, i = 1, 2 and -5 ≤ x k,j ≤ 5, j = 1, 2 respec- tively. The contraction factor is selected to be 0.20. The set shown in front of the Figure 4 is the largest polyhedral contractive set fulfilling these constraints, obtained by the method described in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF]. The dark gray set is the largest ellipsoid and the light gray set at the back is the contractive set obtained by the method described in this paper. The control law for the light gray set is linear and is given as:

u k = -0.0870 -0.0592 -0.0988 0.0653 x k
Table 1 compares the results of the contractive sets obtained by different methods. Dorea-Hennet is the method described in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF] (the set in front in the figure (4)), ellipsoidal is the method which gives the largest ellipsoid fulfilling the constraints (the dark gray set in the figure (4))and VdLF is a variable degree Lyapunov function method developed in this paper. The value of c is 1.22 • 10 -10 .

Example 3

Consider the following linear system, where -5 ≤ x k,j ≤ 5, j = 1, 2, 3 and -2 ≤ u k ≤ 2. The contraction factor is selected to be 0.9. The polyhedral contractive set obtained by the method described in [START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF] is shown in the figure 5 and its volume is 11.9. The largest ellipsoidal set Ω has volume 71.13, and is shown in red in Figure 5. The volume of the set S obtained by the method proposed in the paper, comes out to be 124.36 (green set in figure 5). The set S is the level set of a Lyapunov function of degree 6 and it is obtained using a linear control law.

Discussion and Conclusion

This paper presents a method to obtain a large controlled contractive set of specified degree with the corresponding control law. It was shown that the contractive set can be enlarged by increasing the degree of the contractive set or the degree of the control law. The examples show a significant increase in the size of the contractive set, and hence a corresponding increase in the size of the region where the controller is defined. The resulting controllers are explicit and characterized by low complexity, with a modest memory footprint for the online implementation, and require only very simple computer code.

Currently, the size of SOS programming problems that can be handled is limited, and this will limit both the number of states and the polynomial order that can be handled by the method described in this paper. However, research aiming at increasing the size of SOS problems that can be handled is very active and is making progress, as was documented at the 2017 CDC conference (see, e.g., [START_REF] Zheng | Exploiting Sparsity in the Coefficient Matching Conditions in SOS Programming[END_REF] and [START_REF] Hall | Iterative LP and SOCP-Based Approximations to Semidefinite Programs[END_REF]).

The results in this paper handle nominal linear dynamics only. Conceptually, they are straight forward to generalize to linear systems with polytopic uncertainty, by replicating the design criteria for each extreme point in the polytopic uncertainty set, while keeping the Lyapunov function and the feedback control identical for all extreme dynamics. Simple interpolation can then be used to show that the design criteria also hold for internal points in the uncertainty polytope. However, this will quickly lead to very large SOS problems.

10 :

 10 if solution is feasible & c ≤ c max then 11:

- 5 ≤

 5 x k,j ≤ 5. The specified contraction factor is 0.90. The contractive sets obtained are shown in Figure 1.The set in front is obtained by the Dorea-Hennet procedure[START_REF] Dorea | a,b)-invariant polyhedral sets for linear discrete-time systems[END_REF], it is comprised of 250 hyperplanes. The black set represents the largest ellipsoid obtained with the linear feedback, light and dark gray sets behind black set are obtained using the controller design described in section 3. The light gray set is a level set of the Lyapunov function of degree 8 and the dark gray one is a level set of the Lyapunov function of degree 10. Both of these sets use a linear structure for the control law.

Fig. 1 .

 1 Fig. 1. Contractive Sets of variable degree on V(x) with linear input for Example 1

Fig. 2 .

 2 Fig. 2. State Trajectories for Example 1

Fig. 3 .

 3 Fig. 3. Contractive Sets with varaible input degree for Example 1

Fig. 4 .

 4 Fig. 4. Polyhedral, Ellipsoidal and pre-specified degree Contractive Sets for Example 2

Fig. 5 .

 5 Fig. 5. Controlled contractive sets for example 3.

  Algorithm 1 Algorithm to obtain a larger contractive set Input: A contractive ellipsoid Ω with control law u k . Maximum allowed degree (x deg ) for the Lyapunov function, maximum allowed degree (u deg ) for the control law, and maximum acceptable measure of sub-optimality (c max ). Output: A large contractive set of degree ≤ x deg with control law of degree ≤ u deg . 1: Set j = 0, S j = Ω such that S j = {x k |V (x k ) ≤ 1} and set solution = feasible. LOOP Process 2: while solution is feasible do

	3:

Table 1

 1

	Contractive sets method versus volume
	Method	Degree or Hyperplanes Volume
	Dorea-Hennet	8	75.7677
	Ellipsoidal	2	78.5398
	VdLF	8	87.2972

For subsequent developments to make sense, the controller Kuc should clearly be designed such that the unconstrained closed loop system is γ-contractive, i.e., such that max |eig(A + BK)| ≤ √ γ.

The work leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement no 607957 (TEMPO).