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.