Self-tuning model predictive control for wake flows

Abstract This study presents a noise-robust closed-loop control strategy for wake flows employing model predictive control. The proposed control framework involves the autonomous offline selection of hyperparameters, eliminating the need for user interaction. To this purpose, Bayesian optimization maximizes the control performance, adapting to external disturbances, plant model inaccuracies and actuation constraints. The noise robustness of the control is achieved through sensor data smoothing based on local polynomial regression. The plant model can be identified through either theoretical formulation or using existing data-driven techniques. In this work we leverage the latter approach, which requires minimal user intervention. The self-tuned control strategy is applied to the control of the wake of the fluidic pinball, with the plant model based solely on aerodynamic force measurements. The closed-loop actuation results in two distinct control mechanisms: boat tailing for drag reduction and stagnation point control for lift stabilization. The control strategy proves to be highly effective even in realistic noise scenarios, despite relying on a plant model based on a reduced number of sensors.


Introduction
Prediction and control of fluid flows to pursue a specific objective is a highly compelling research area (Gad-el Hak 2000).Flow control offers wide-ranging practical applications in diverse fields, including vehicle dynamics, aircraft, and marine transportation, meteorology, energy production from water and wind, combustion, and chemical processes, and more (Duriez et al. 2017).The goals of fluid flow control generally encompass, among others, drag reduction, control of separation and transition, lift or mixing enhancement.In recent years, drag reduction has received considerable attention due to its notable impact on the environmental footprint of transportation means (Green 2003).
In the last decades, active flow control has garnered increasing attention.This technique can be implemented in open-loop and closed-loop configurations.The former involves predetermining the actuation law irrespective of the system state, thus simplifying its application.Notable examples include wake control of bluff bodies (Blackburn & Henderson 1999;Cetiner & Rockwell 2001;Thiria et al. 2006;Parkin et al. 2014) and open-cavity flows (Sipp 2012;Little et al. 2007;Nagarajan et al. 2018), among others.However, open-loop control's effectiveness is limited in unstable flow stabilization, responding to changes in the controlled system parameters, or dealing with external disturbances.On the contrary, closed-loop implementations (also referred to as reactive control) involve feeding the control law by the knowledge of the state.This approach offers greater flexibility and adaptability (Brunton & Noack 2015).Experimental evidence demonstrates the superior performance of closed-loop over open-loop control, see e.g.Pinier et al. (2007) or Shimomura et al. (2020).
The identification of control laws requires adequate knowledge of the system dynamics and † Email address for correspondence: lmarra@pa.uc3m.esarXiv:2401.10826v1[physics.flu-dyn]19 Jan 2024 its response to control inputs.In fluid dynamics, model-based techniques have traditionally been utilized to obtain this information, proving successful in various scenarios (Kim & Bewley 2007).Examples of applications include transition delay in spatially evolving wallbounded flows (Monokrousos et al. 2008;Chevalier et al. 2007;Tol et al. 2019), cavity flow control (Rowley & Williams 2006;Illingworth et al. 2011), separation control on a low Reynolds number airfoil (Ahuja et al. 2007), wake stabilization of cylinders (Schumm et al. 1994;Gerhard et al. 2003;Tadmor et al. 2011), skin-friction drag reduction (Cortelezzi & Speyer 1998;Lee et al. 2001;Kim 2011) and heat transfer control (Castellanos et al. 2022b).However, the identification of efficient analytical control laws faces an important challenge in the presence of complex nonlinear multiscale dynamics.
In recent years, model-free techniques have gained popularity, driven by advancements in hardware and the increasing efficiency of data-driven and machine-learning (ML) algorithms.Examples of model-free techniques include genetic algorithms in jet mixing optimization (Koumoutsakos et al. 2001;Wu et al. 2018), wake flows (Poncet et al. 2005;Raibaudo et al. 2020), separation control (Gautier et al. 2015), and combustion noise (Buche et al. 2002).Reinforcement learning (RL) has also recently gained popularity, with successful applications in the control of bluff body wakes (Rabault et al. 2019;Fan et al. 2020;Castellanos et al. 2022a) and natural convection (Beintema et al. 2020).Despite the encouraging results of such model-free techniques, their effectiveness is limited by the need for large datasets.
Within model-based techniques, model predictive control (MPC) offers interesting features to deal with the challenges of fluid flow control.MPC is based on the idea of receding horizon control.It found application in the industry since the 1980s (Qin & Badgwell 2003), in particular with extensive use in refineries and the petrochemical industry (Lee 2011).MPC has demonstrated excellent performance in controlling complex systems with constraints, strong non-linearities, and time delays (Henson 1998;Allgöwer et al. 2004;Camacho & Alba 2013;Grüne & Pannek 2017).Therefore, it is particularly appropriate for complex systems that challenge traditional linear controllers (Corona & De Schutter 2008).The method requires the identification of a model of the system dynamics capable of predicting its behavior under exogenous inputs.The optimal control is determined through the iterative solution of an optimization problem within a prediction window, aiming to minimize a userdefined cost function that considers the distance of the system state from the control target.Moreover, MPC allows for the straightforward implementation of hard constraints, such as hardware limitations, distinguishing it from classical control approaches.MPC has been successfully applied in the control of complex fluid systems, see e.g. the works of Collis et al. (2000); Bewley et al. (2001); Bieker et al. (2020); Sasaki & Tsubakino (2020); Morton et al. (2018); Peitz et al. (2020).A crucial aspect of MPC implementation is achieving a proper balance among the terms of the loss function.The user needs to select weights (referred to as hyperparameters) for the loss, considering factors like closeness to the target, cost of the action, and other application-tailored constraints.This choice has a clear impact on the final performance.In flow control applications, this process traditionally relies on user experience, which poses the risk of suboptimal choices.
Bayesian optimization (BO) or RL techniques have demonstrated excellent results in hyperparameter tuning, particularly in the fields of autonomous driving and robotics (Edwards et al. 2021;Fröhlich et al. 2022;Bøhn et al. 2023).A comprehensive review in this area can be found in Hewing et al. (2020).However, in the application of nonlinear MPC to flow control, examples are scarce, and the choice of MPC parameters is often guided by trial-error and intuition.This approach risks falling into suboptimal configurations that may not adequately account for the different degrees of fidelity in the terms involved in the loss function.This issue is particularly relevant in fluid mechanics, where the uncertainty in predicting the plant behavior and the measured state/control actions should play a role in the parameter selection process.Unfortunately, an analytical formulation is elusive in most cases.
Moreover, as a closed-loop strategy, the implementation of MPC requires feedback, consisting of time-sampled measurements of a feature of the system to be controlled.In real control scenarios, this sampling is often affected by measurement noise, which can compromise control decision-making.Thus, suitable smoothing techniques are necessary to enhance noise robustness.In time series analysis, a nonparametric statistical technique called local polynomial regression (LPR) proves particularly effective in this task.LPR estimates the regression function of sensor outputs and their time derivatives without assuming any prior information.Applications of LPR for control purposes are described in works such as Steffen et al. (2010) or Ouyang et al. (2018).
In this article, we propose a fully automatic architecture that self-tunes control and optimization process parameters with minimal user input.Our MPC framework adapts to different levels of noise and/or limited state knowledge.The methodology builds upon offline black-box optimization via Bayesian methods for hyperparameter tuning.Furthermore, we discuss the robustness enhancement to noise using an online local polynomial regression.The effectiveness of the control algorithm is evaluated through its application to the control of the wake of the fluidic pinball (Deng et al. 2020) in the chaotic regime.Although not strictly needed, plant identification is also data-driven.In this work, nonlinear system identification is performed using the sparse identification of nonlinear dynamics with control (SINDYc, Brunton et al. 2016b).
The paper is organized as follows.Section 2 provides a description of the methodology, emphasizing the mathematical tools and the MPC framework employed.Additionally, this section includes specific details regarding the chosen test case for illustration purposes.The results of the control application, along with their interpretation are provided in § 3. Finally, the conclusions are discussed in § 4.

Methodology
This section presents the backbone of the MPC-based control algorithm, with a detailed description of the mathematical tools involved in it.Figure 1 includes a diagram illustrating all the required steps for its implementation.In addition, algorithm 1 is introduced to give more detail on the procedure.
The main block in the schematics represents the MPC algorithm, following the approach proposed by Kaiser et al. (2018).The main novelty in this module is the robustness enhancement by online filtering with LPR.This is particularly useful when the plant dynamics is modeled with time-delay coordinates or their derivatives.LPR is directly applied to past sensor data for online implementation.
The roadmap suggests several necessary steps before implementing the control.First, a training dataset is generated.The dataset consists of the time series of the state dynamics under different exogenous inputs.This data can be collected using various methods, including simulations or experiments.The system state and exogenous inputs should be defined based on the specific system being controlled.Second, a plant model is defined to predict the system behavior.In this work, we use a data-driven nonlinear system identification.The final step, before the implementation of the control, focuses on tuning the parameters that define the MPC cost function.Control performance is significantly influenced by their selection.The self-tuning of the hyperparameters is the core of our framework.This tuning is carried out using a Bayesian optimization algorithm.
An important aspect of the proposed method is the need for minimal user interaction.Indeed, the two main inputs are the reference setpoint (i.e. the control target), and a cost

Nonlinear system identification
Reduced-force model with relatively simple dynamics (Deng et al. 2020(Deng et al. , 2022)), which evolve on a low-dimensional attractor.The extension of the method to the control of turbulent flow cases is left for future investigations, where dealing with measurement noise in control sensors and selecting a reduced model of the system is more complex.

Control
Another aspect that requires further attention and improvement compared to the proposed method is the selection of the predictive model.In future work, it is suggested to test the prediction performance using various other methods, including cluster-based methods and NARX, among others.
Funding.The authors acknowledge the support from the research project PREDATOR-CM-UC3M.This project has been funded by the call "Estímulo a la Investigación de Jóvenes Doctores/as" within the frame of the Convenio Plurianual CM-UC3M and the V PRICIT (V Regional Plan for Scientific Research and Technological Innovation).function for the MPC algorithm.The weight of the different contributions in the cost function will be determined in the MPC-tuning.BO automatically adjusts to different levels of noise and uncertainties.This greatly enhances the usability and adaptability of the framework to different systems.

Model predictive control
This section describes the MPC implementation.It is assumed to have a time-evolving process whose complete description relies on   ≥ 1 parameters.These are included in a state vector, denoted at a given instant  as  ≡ () =  1 (), . . .,    ()

′
. The evolution in time of the process can be influenced by the choice of   ≥ 1 exogenous parameters.These are included in an input vector  ≡ () =  1 (), . . .,  dynamics is described by the following set of equations: where (  ) =   and  is the function characterizing the system's temporal evolution.The process is considered as controlled.This means that for each time , the input vector () can be selected in order to manipulate the system dynamics according to a specific objective.More specifically, the aim is to control   ≥ 1 features of the dynamical system, included in the vector  ≡ () =  1 (), . . .,    () ′ .Note that the vector  is dependent on the system state.In this context, it is assumed that the target features are part of the state vector itself, although this may not always be applicable.The objective is to achieve control over the vector  by ensuring that it closely tends to a desired reference  * ∈ R   over time.MPC can be used for setpoint stabilization, trajectory tracking, or path following (see Raković & Levine 2018, pp. 169-198), depending on the choice of  * .
For the purpose of a control application, ,  and  are sampled over a discrete time vector, equispaced with a fixed time interval   .This discrete-time representation is essential to determine how often the exogenous input is updated.The input is assumed to be constant between consecutive timesteps of the control.The implementation of the MPC follows a series of sequential procedures, as can be seen in algorithm 1.An illustration of the process is provided in figure 2.
Firstly, the control process starts from a time instant   , where a measurement   of the state vector is available.It is assumed that the entire vector of target features is observed.A conditional prediction of the state vector is obtained by a model of the dynamics.This prediction under a given input sequence, referred to as â + |  , is generated within a prediction window  + ,  = 1, . . .,   .Consequently, a prediction of the target features vector ĉ + |  is obtained.
The optimal input sequence {    }   =1 can be determined in a control window  + ,  = 1, . . .,   , by minimizing a cost function J   : R   → R + .A common choice for the cost function is where Q ∈ R   ×   and R  , R Δ ∈ R   ×   are positive and semi-positive definite weight matrices, respectively.In addition, ∥ ∥ 2 M =  ′ M  represents the weighted norm of a generic vector  with respect to a symmetric and positive definite matrix M, where  ′ denotes the transposition of the vector .The variable Δ + |  , denotes the input variability in time, that is  + |  −  +−1|  .In the definition of the cost function, the errors in state predictions with respect to the reference trajectory, i.e. ĉ + |  −  * , are penalized, as well as the actuation cost and variability.In this paper, the aforementioned weight matrices are assumed to be diagonal, thus: (2.3) Once the optimization problem in (2.2) is solved, only the first component of the optimal control sequence,  +1 =    1 , is applied.The optimization is then reinitialized and repeated at each subsequent timestep of the control.Generally, the prediction window covers a wider range than the control window (  ≥   ).The control vector is considered constant beyond the end of the control window, as discussed in Kaiser et al. (2018).
Hard constraints on the input vector can readily be incorporated into the optimization process.At each timestep, the optimal problem must guarantee that  ∈ B, where B is generally determined by the control hardware.The set of allowable inputs is as follows: A critical issue of this control technique concerns the choice of parameters.Many of these can be selected based on the physics of the system to be controlled or appropriate control hardware limits, such as the actuation constraints or the control timestep.The selection of parameters involved in equation (2.3), as well as the length of the prediction/control window, is traditionally tuned by trial and error.In § 2.1.4we introduce our hyperparameter self-tuning procedure.

Nonlinear system identification
The optimization loop of MPC requires an accurate plant model for predicting the system's dynamics.The choice of the predictive model typically involves a trade-off between accuracy and computational complexity.In this work, we adopt a data-driven sparsity-promoting technique.The framework was illustrated in the work by Brunton et al. (2016a) and later applied in MPC by Kaiser et al. (2018) for a variety of nonlinear dynamical systems.It must be remarked that the self-tuning framework proposed here can easily be adapted to accommodate a different plant model, either analytical or data-driven, including deep-learning models.
SINDy is a method that identifies a system of ordinary differential equations describing the dynamical system.In particular, the version described in this work corresponds to the extension of the SINDy model with exogenous input (SINDYc, Brunton et al. 2016b).Considering a system of ordinary differential equations such as the one shown in (2.1), SINDYc derives an analytical expression of  from data.This process requires a dataset comprising the time series of the state vector  and the exogenous input .This method is based on the idea that most physical systems can be characterized by only a few relevant terms, resulting in governing equations that are sparse in a high-dimensional nonlinear function space.The resulting sparse model identification aims to find a balance between model complexity and accuracy, preventing overfitting of the model to the data.
To derive an expression of the function  from data, a discrete sampling is performed on a time vector, yielding  snapshots at time instances   ,  = 1, . . .,  for the state vector   = (  ), its time derivative   = (  ), and the input signal   = (  ).The data obtained are arranged into three matrices (2.4) A library of functions  is then set as: where 1  is a column vector of  ones, A • denotes the   ℎ row of the matrix A and H ⊗ K denotes the Kronecker products of H and K .The choice of functions to be included in the library is typically made by the user.This procedure is often guided by experience and prior knowledge about the dynamics of the system to be controlled.This issue introduces some limitations that are later discussed in §4.
The system data can be then assumed to be generated from the following model: =  1 , . . .,    represents the matrix whose rows are vectors of coefficients determining which terms of the right-hand side are active in the dynamics of the   ℎ state vector component.This matrix is sparse for many of the dynamical systems considered.It can be obtained by solving the optimization problem: where  is called sparsity-promoting coefficient, A • denotes the   ℎ column of the matrix A. ∥•∥ 1 and ∥•∥ 2 are the  1 and  2 norms, respectively.It must be remarked that  must be optimized to reach a compromise between parsimony and accuracy.Nonetheless, since this block of the MPC framework can easily be replaced by other strategies,  will not be included in the self-tuning optimization illustrated in § 2.1.4.Finally, the system can be described by: where in this case (, ) is a vector that takes into account the same functions included in the library in (2.5).The SINDYc-based model is utilized to make predictions of the dynamics involved in the control process starting from a specific state's initial condition.The steps described in this section to obtain the plant model are also condensed in the nonlinear system identification step in algorithm 1.

Local Polynomial Regression
Effectively estimating system dynamics is crucial for various control techniques, particularly in MPC, where sensor feedback is utilized to make informed decisions.The accuracy of system dynamics estimation becomes paramount in MPC as it relies on sensor information to predict the behavior of the controlled system.However, the presence of noise can hinder accurate estimation, necessitating the implementation of noise mitigation techniques.This challenge falls within the broader context of time series analysis, as sensor outputs represent discrete samples over time of a process variable.To address measurement noise, LPR emerges as a robust, accurate, and cost-effective nonparametric smoothing technique.Applying LPR to sensor output data enhances the reliability and accuracy of estimating the current state of the system under varying noise conditions.In this section, we provide a succinct overview of the formulation of LPR.The notation convention involves denoting random variables with capital letters, while deterministic values are represented in lowercase.
The objective is to predict or explain the response  (sensor output) using the predictor  (time) from a sample {(  ,   )}  =1 using the following regression model: where  is the regression function, {  }  =1 , are zero mean and unit variance random variables and  2 is the point variance.For simplicity, it is assumed that the model is homoscedastic, i.e. that the variance is constant.From a practical point of view,  also represents the measurement of noise intensity.
The reason why LPR is used in this context is that it allows for a fast joint estimation of the regression function and its derivatives without making any prior.By making only a few regularity assumptions, it is characterized by a good flexibility that allows to model complex relations beyond a parametric form.The working principle of the LPR is to approximate the regression function locally by a polynomial of order , performing a weighted least squares regression using data only around the point of interest.
It is assumed that the regression functions admit derivatives up to the order  + 1 and thus that it can be expanded in Taylor's series.For close time instants  and   the unknown regression function can be approximated by a polynomial of order , then The vector of parameters  = ( 0 , . . .,   ) ′ can be estimated by solving the minimization problem: where ℎ is the bandwidth determining the size of the local neighborhood (also called smoothing parameter) and  ℎ () = 1 ℎ (  ℎ ) with  a kernel function assigning the weights to each observation.Once the vector β has been obtained, an estimator of the   ℎ derivative of the regression function is m() () = !β .The solution of the optimization in (2.9) is simplified when the methodology is presented in matrix form.To that end, T ∈ R ×( +1) is the design matrix while the vector of responses is  = ( 1 , . . .,   ) ′ and W = diag( ℎ ( 1 − ), . . .,  ℎ (  − )).
According to this notation and assuming the invertibility of T ′ WT , the weighted least squares solution of the minimization problem expressed in (2.9) is and thus the estimator of the local polynomial is: being   ∈ R ( +1)×1 a vector having 1 in the   ℎ entry and zero elsewhere and The critical element that determines the degree of smoothing in LPR is the size of the local neighborhood.The selection of ℎ determines the accuracy of the estimation of the regression function: too large values lead to a large bias in the regression, while too small values increase the variance.Its choice is generically the result of a trade-off between the variance and the estimation bias of the regression function.There are several methods to select the bandwidth to be used for LPR.One possible approach is to choose between a global bandwidth, which is common to the entire domain and is optimal for the entire range of data, or a local bandwidth that depends on the covariate and is therefore optimal at each point of the regression function estimation.The latter would allow for more flexibility in estimating inhomogeneous regression functions.However, in the proposed framework, a global bandwidth was chosen due to its simplicity and, more importantly, its reduced computational cost.In the presented approach, the global bandwidth is chosen using the leaveone-out-cross-validation (LOOCV) methodology.Specifically, this parameter corresponds to the one that minimizes the following expression: where mℎ,−  (  ) denotes the estimation of the regression function while excluding the   ℎ term.For further information regarding the optimal bandwidth selection see the works of Wand & Jones (1994) and Fan & Gijbels (1996).Additionally, there exist other methodologies that utilize machine learning techniques to select the local bandwidth such as in Giordano & Parrella (2008) where neural networks are used.
Once the smoothing parameter has been set, it remains to choose the weighting function and the degree of the local polynomial, although these two have a minor influence on the performance of the LPR estimation.A common choice for the first is the Epanechnikov kernel.As for the degree of the local polynomial, there is a general pattern of increasing variability according to which, to estimate  () (), the lowest odd order is recommended, i.e.  =  + 1 or occasionally  =  + 3 (Fan & Gijbels 1996).
It is also worth noting that due to asymmetric estimation within the control algorithm, boundary effects arise, leading to a bias in the estimation at the edge.These effects, discussed in more detail in Fan & Gijbels (1996), disappear when using local linear regression ( = 1).

Hyperparameter automatic tuning with Bayesian optimization
MPC relies on a specific set of hyperparameters to precisely define the cost function in (2.2) for the selection of the optimal action over the control window.A new functional, based on the global performance of the control, is defined.Bayesian optimization is employed to find the hyperparameter vector that maximizes the control performance.By adopting this approach, a more efficient and effective MPC-based control framework may be achieved.The hyperparameters are indeed adapted to different conditions of measurement noise and/or model uncertainty by the BO process.This proposed framework is practically independent of the user.
The parameters under consideration include i) the components of the weight matrices, as presented in (2.3), which penalize errors of the state vector with respect to the control target trajectories, input cost, and input time variability and ii) the length of the control/prediction windows, assumed equal for simplicity.All aforementioned control parameters are included in a single vector denoted as  ∈ R   , where   = 2  +   + 1.A further reduction in the number of control parameters can be also considered by imposing that the components of R  and R Δ related to the rear cylinders of the fluidic pinball are identical under a flow symmetry argument, as in Bieker et al. (2020).Therefore, in this case it is stated that   2 =   3 =   2,3 and  Δ 2 =  Δ 3 =  Δ 2,3 and the total number of parameters is reduced to   = 2(  − 1) +   + 1. Section 3 also provides several results that justify this choice.
Note that control results depend on the selection of the vector .Consequently, an offline optimization process can be implemented to select the optimal value of , which maximizes the control performance.Consider using a specific realization of the hyperparameter vector  to control the system.The state vector measure is indirectly dependent on the choice of hyperparameter vector and is denoted here by s().By running the control on a discrete-time vector of   timesteps, the cost function can be defined as follows: (2.10) Note that the user selects the target to be optimized, specified by the cost function J  .In addition, in absence of measurement noise, s() is the ideal measure of the target feature vector.Otherwise, the LPR technique is applied to the state vector measure.In this study, the parameters are optimized to minimize the quadratic error between the controlled variables and the user-defined target.Thus, the optimal hyperparameter vector can be obtained by solving the problem where  ⊂ R   is the search domain.Solving the optimization in (2.11) is computationally costly due to the expensive sampling of the black-box cost function J  .Each sample of J  is obtained via application of the MPC-based control over   timesteps.In this framework, BO serves as an efficient method to address this problem, offering an algorithm to search the minimum of the function with a high guarantee of avoiding local minima and characterized by rapid convergence.Indeed, BO has shown to be an efficient strategy particularly when   ≤ 20 and the search domain  is a hyper-rectangle, that is where   and   are the lower and upper bound vectors of the hyper-rectangle, respectively.
In order to find the minimum of the unknown objective function, BO employs an iterative approach, as shown in the MPC-tuning step of algorithm 1.It utilizes a probabilistic model, typically a Gaussian process (GP), to estimate the behavior of the objective function.At each iteration, the probabilistic model incorporates available data points to make predictions about the function behavior at unexplored points in the search space.Simultaneously, an acquisition process guides the samplings by suggesting the optimal locations in order to discover the minimum point.A specific function (called acquisition function) is set to balance exploration, by directing attention to less-explored areas of the search domain, and exploitation, by concentrating on regions near potential minimum points.Finally, the BO iterations continue until a stopping criterion is met, such as reaching a maximum number of iterations or achieving convergence in the search for the minimum.Thus, denoting as () the acquisition function selected, the point where to sample next in the iterative approach, that is  + , can be obtained by solving The expected improvement (EI) is used as an acquisition function (Snoek et al. 2012), which evaluates the potential improvement over the current best solution.
The subsequent discussion will center on the probabilistic model utilized in Bayesian optimization.Specifically, a prior distribution, which corresponds to a multivariate Gaussian distribution, is employed.Consider the situation where problem (2.11) is tackled using Bayesian optimization.A Gaussian process regression is required for the functional J  based on the available observations of the function up to the given iteration.
Since J  is a Gaussian process, for any collection of  points, included in  =  1 , . . .,   ′ , then the vector of function samplings at these points, denoted as  = J  ( 1 ), . . ., J  (  ) ′ is multivariate-Gaussian-distributed, then: being  = ( 1 ), . . ., (  ) ′ the mean vector and K the covariance matrix whose    ℎ component is    = (  ,   ), with  a mean function and  a positive definite kernel function.For simplicity, the mean function is assumed to be null.
In order to make a prediction of the value of the unknown function at a new point of interest  * , denoted as J  * , conditioned on the values of the function already observed, and included in the vector , the joint multivariate Gaussian distribution can be considered: where  * * = ( * ,  * ), and the vector  * = ( 1 ,  * ), . . ., (  ,  * )

′
. This equation describes how the samples  at the locations  correlate with the sample of interest J  * , whose conditional distribution is: . Equation (2.12) provides the posterior distribution of the unknown function at the new point where the sample has to be performed.It then furnishes the surrogate model used to describe the function J  in the domain for the search of the minimum point.

The fluidic pinball
The proposed framework is tested on the control of the two-dimensional viscous incompressible flow around a three-cylinder configuration, commonly referred to as fluidic pinball  The dynamics of the wake past the fluidic pinball is obtained through a two-dimensional direct numerical simulation (DNS) with the code developed by Noack & Morzyński (2017).The flow is described in a Cartesian reference system in which the  and  axes are in the streamwise and crosswise directions, respectively.The center of the Cartesian reference system coincides with the mid-point of the rightmost bottom and top cylinder and the computational domain, that is Ω, is bounded in (−5, 20) × (−5, 5).The position in the reference system is then described by the vector  = (, ) =  1 +  2 , where  1 and  2 are respectively the unit vectors in the directions of the  and  axes and the velocity vector is assumed to be  = (, ).The constant density is denoted by , the kinematic viscosity of the fluid by .All quantities used in the discussion are assumed non-dimensionalized with cylinder diameter, freestream velocity, and fluid density.The Reynolds number is defined as Re  =  ∞   .A value of Re  = 150 is adopted, which is sufficiently large to ensure a chaotic behavior, although still laminar.The 2D solver has already been used in previous work at this same Re  (Wang et al. 2023).The boundary conditions comprise a far field condition ( =  ∞  1 ) in the upper and lower edges, a stress-free one in the outflow edge, and a no slip-condition on the cylinder walls, which in the absence of forcing becomes  = 0.
The DNS allows forcing by independent rotation of the cylinders.To this purpose, a velocity with module |  | is imposed at the cylinder surface.Positive values of   are associated with counterclockwise rotations of the cylinders.In the remainder of the paper, a reference timescale is set as the convective unit (..), i.e. the timescale based on the freestream velocity and the cylinder diameter.The lift coefficient is defined as ∞  , being   the total lift force applied to the cylinders in the direction of -axis and the same quantity is applied in the scaling of   , force applied to the cylinders in direction of -axis, to obtain the total drag coefficient   .The DNS uses a grid of 8633 vertices and 4225 triangles accounting for both accuracy and computational speed.A preliminary grid convergence study at Re  = 150 identified this as sufficient resolution for errors of up to 3% in the free case and ≈ 4% in the actuated case in terms of drag and lift.

The control approach
The aim of the control is to achieve a reduction in the overall drag coefficient of the fluidic pinball, while also controlling the lift coefficient so that the latter has reduced oscillations with a zero average value.The vector of target features is therefore composed only of the total lift and drag coefficients,  = (  ,   ) ′ .The target vector is thus composed of null components  * = (0, 0) ′ .Since the cost function in (2.2) is quadratic, a null target vector will penalize both mean values of   and   and their oscillations.
To this purpose, the tangential velocities of the three cylinders are tuned respecting the implementation limits, here chosen equal to   ∈ [−1, 1] and Δ  ∈ [−4, 4] , ∀ = 1, 2, 3.The model plant has a state vector composed by   and   and their time derivatives, respectively   and   (i.e. = (  ,   ,   ,   ) ′ ).A state vector based solely on global variables does not require adding intrusive probes for flow estimation.Similar state vector choices that incorporate aerodynamic force coefficients and temporal derivatives are made in works Nair et al. (2019) and Loiseau et al. (2018).While this approach is often effective in separated flows, it might be challenged at higher Re flows and of unfeasible application in other flow configurations.
Feedback to the control involves measurement of drag and lift forces, so at each time instant   , it is considered   =   .This approach also facilitates a reduction in the complexity of the predictive model, thereby speeding up the control process.Indeed, a compact state vector is particularly desirable to reduce the computational cost of the iterative optimization problem over receding horizons of MPC.
In the present work, noise in lift and drag measurement is also considered.Under this assumption, the model in equation (2.8) is applied.The noise intensity  is assumed to be constant over time and given as a percentage of the full-scale measured drag and lift coefficients in conditions without actuation.Therefore in the case of measurement noise in the sensors, the response variable to which the LPR is applied corresponds to the measurements of the drag and lift coefficients of the fluidic pinball.The LPR enables us to obtain estimates of their regression function in output from the sensors but also of their time derivatives, allowing us to use this information as control feedback.
The MPC-optimization problem is solved every control timestep (so every   ) to update the exogenous control input.During the time between two consecutive samples, the control input to the system remains constant.Thus, the sampling timestep should be chosen small enough to ensure a good closed-loop performance, but not too small to avoid an excessive computational cost.In this work, it is set at   = 0.5 .., i.e. sufficiently small compared to the shedding period of the fluidic pinball wake (denoted as  ℎ ) and not too high to affect control performance.The reason why   was not included in the control tuning concerns the difficulty of defining the search domain for realistic applications.Imposing bounds on this parameter requires an estimate of the computational cost of the MPC-based control, and this procedure is postponed to future experimental applications.The method used to optimize the control action in the MPC framework is sequential quadratic programming with constraints.The optimization was carried out with a built-in function of MATLAB.The stop criteria are set at a maximum number of iterations of 500 and a step tolerance of 1−6.
At each measurement taken during the control process, in the presence of noise in the   and   sensors, the LPR technique is applied to a time series of outputs of length corresponding to the characteristic shedding period of the fluidic pinball wake, as observed in § 3. LPR acts asymmetrically, i.e.only past output sensor data is available.
Parameter tuning is performed by evaluating the cost function J  over a time history of the controlled variables of one shedding period of the fluidic pinball.The transient phase is excluded.Finally, as a stopping criterion for the search for the optimal hyperparameter vector, a maximum number of 100 iterations in Bayesian optimization was set.This proved sufficient to reach convergence in the tuning process, as shown in § 3.2.

Results
In this section, the results of the current work are presented.Subsection 3.1 outlines the prediction performance of the SINDYc-based force model.Furthermore, the prediction performance enhancement when using LPR in the presence of measurement noise is illustrated.In § 3.2 the effectiveness of hyperparameters tuning in the different control scenarios is examined.Lastly, in § 3.3, the outcomes of applying the MPC-based control algorithm to the fluidic pinball wake for drag reduction purposes under ideal measurement conditions and in the presence of sensor noise are presented.An additional discussion is provided to justify the choice of symmetry in the parameters of the rear cylinders of the fluidic pinball.

Predictive model
This subsection presents the performance evaluation of the predictive model employed in the control framework.SINDYc is applied on a training dataset consisting of a time series of the system state, along with predetermined actuation laws.The chosen open-loop laws are composed of step signals for the three-cylinder velocities   , encompassing all possible combinations of the velocities of the three cylinders within control constraints.Specifically, velocities ranging from −1 to 1 with an increment of 0.5 are explored, resulting in a total of 5 3 combinations.For each explored velocity combination, the actuation steps are long enough (≈ 55 ..) to allow reaching the respective steady state.The total length of the generated dataset is 6950 ...The actuation time series is subsequently smoothed to include the transient dynamics in the force model.It must be noted that this approach is effective for this test case whose dynamics evolve on a clearly defined low-dimensional attractor.Openloop training for plant identification has been shown effective in separated flows (see e.g.Kaiser et al. 2017;Bieker et al. 2020).More complex nonlinear dynamical systems might require different strategies for adequate training of the modeling plant.
Additionally, a library of polynomial functions up to the second order is chosen for the sparse regression with SINDYc.
Figure 4 shows the prediction performance of the plant model.The presented results were generated by performing predictions on a testing dataset with preassigned smooth law of the rotational speeds of the three cylinders.Predictions of 4 ..using initial conditions at random points of the validation dataset, and for a statistically significant number of times were performed.The resulting time series of each prediction were then used to calculate their respective errors ( ê  for   and ê  for   ).The plot illustrates the trend of the average and the confidence region for  of the prediction errors, respectively normalized with respect to the standard deviation of   and   in the free response, here denoted as (  0 ) and (  0 ), respectively.
The panels on the left show the prediction performance calculated from an initial condition obtained in the absence of measurement noise in the sensors.The average value of the normalized error distribution maintains values close to zero for any length of the prediction window analyzed.Furthermore, for short prediction lengths (up to ≈ 3 ..) the confidence bounds for  are still confined within a narrow region.In this context, as they are not directly measurable, the time derivatives of the aerodynamic coefficients are calculated using finite differences.This approach, however, is highly sensitive to measurement noise.The nonparametric smoothing technique LPR allows joint estimation of the output data from the sensors and their respective time derivatives.The panels on the right show the benefits brought by initial-condition setting with LPR to the prediction performance in the presence of a sensor measurement noise of 1%.These plots show a comparison of the statistical error distributions for two distinct cases: in one the initial condition corresponds to the LPR estimation of the data with noise, in the other, it is not used.It is of considerable interest to observe that The prediction with LPR performs similarly to the case without noise.The figure above discussed could be used to choose the MPC prediction window length as it provides an idea of the uncertainty propagation in the prediction horizon.However, since this parameter has a significant effect on control results in terms of performance and computational cost, it is left to automatic selection through Bayesian optimization, so as to also take into account the effects of measurement noise in its selection.In addition, figure 5 provides further insights into the advantages of LPR.In this case the estimation of   and   is tested on a case with external forcing and under increasing noise level.The trend in the LPR estimation error of both   and   is illustrated as the noise increases.The estimation errors are evaluated using a root mean square error (RMSE) of the estimate compared to the ideal data without noise.
The results demonstrate the effectiveness of LPR in estimating   and   with high accuracy even in the presence of high levels of noise.Thus, the benefits brought by LPR enable the control feasibility even in high-noise scenarios.

Bayesian optimization for tuning control parameters
This subsection presents the results of the hyperparameters tuning using the Bayesian optimization algorithm.Various control scenarios are analyzed, including clean, low, and high noise level cases.
The set of hyperparameters to be optimized is composed of the elements of the weight matrices in (2.3) and the length of the control/prediction window, used in the definition of the cost functional in (2.2).Specifically, the components of the matrix Q are optimized within the interval [0.1, 5], the terms of R  and R Δ in [0.1, 10], and   is optimized within [1,4].
Figure 6 displays the trend of the minimum observed J  samplings as the tuning process iterations progress for some of the analyzed control cases.The results indicate that the optimization function (J  ) reaches convergence in less than 30 iterations for almost all cases.The tuned control parameters of each case are presented in table 1.The table also provides the average values of   and   ( C and C ), as well as their corresponding standard  deviations ((  ) and (  )), during the control phase.These values are computed over a shedding time interval.
It must be remarked that the final set of hyperparameters for each case is run-dependent.This is mostly to be ascribed to the weak dependence of J  to some of the imposed constraints in the MPC loss function.Nonetheless, the differences in terms of drag reduction have been observed to be smaller than ±1.1% in different runs.
Figure 7 presents a plot of the cost function J  and its contributions due to the control of   (denoted as J 1  ) and   (denoted as J 2  ), varying  1 and  2 , and then   1 and   2 =   3 =   2,3 under 1% measurement noise.The two contributions of J  can be obtained by considering the cost functional in equation ( 2 provided for all the analyzed simulation cases, with a particular focus on the tuning effects due to the intensity of the sensor measurement noise, the use of the LPR technique in the sensor outputs, and the imposition of symmetry in the control parameters.The asterisk ( * ) denotes control cases in which the tuning is performed on all control parameters, without imposing the symmetry condition on the parameters of the rear cylinders of the fluidic pinball.The table also provides the case with all parameters equal to 1 (except for the prediction/control window set at 3..) the component for  = 1 (to obtain J 1  ) and for  = 2 to obtain J 2  .The presented cost function plot is useful for interpreting the convergence parameters of the optimization process according to BO.The term J 2  has a minor influence on the optimization process, being approximately two orders of magnitude smaller than J 1  .Therefore, the optimization process is mainly driven by the influence of the control parameters on reducing the   of the fluidic pinball.
As  1 increases, both J  and J 1  functions decrease sharply at first and then reach a plateau.Conversely, J 2  shows the opposite behavior.A good   control would be achieved with a low  1 and high  2 , but this would lead to a poor drag reduction.On the other hand, the choice of  2 has little impact on control performance, as the cost function J  remains nearly constant with its variation.
By observing table 1 and considering the behavior of the cost function, it is justified why almost all proposed control cases have  1 very close to the maximum achievable value, while there is greater variability in the selection of  2 .
Similarly, due to the fact that the drag reduction drives the optimization process, the parameter that has the greatest influence on J  is   2,3 .A lower value of   2,3 would assign less weight to the actuation cost of the rear cylinders of the fluidic pinball, allowing for greater rotation intensity.As it will be explained in § 3.3.1, the main contributors to drag reduction are indeed the rear cylinders.On the other hand, the parameter   1 has little impact on the overall control performance since the front cylinder is responsible for the stabilization of the lift coefficient.These latter considerations justify that the optimum   2,3 is achieved at low values, close to the lower possible limit, while greater variability is observed in the optimization of the   1 parameter.

Application of the MPC-based control algorithm
This subsection shows the results of applying the control algorithm to the wake of the fluidic pinball.It is specified that all results reported from now on are obtained by imposing symmetry in the parameters of the rear cylinders.).Contour plots are presented as functions of  1 and  2 (top graphs) and   1 and   2 =   3 =   2,3 (bottom graphs).The control case considers a 1% measurement noise and symmetric weighting coefficients for the rear cylinders of the fluidic pinball.In the graphs, the remaining components of the hyperparameter vector  are fixed at their optimal values.

Control simulations without measurement noise
Before discussing the effects of MPC-based control, a brief description of the behavior of the unforced wake of the fluidic pinball in the laminar chaotic regime, characteristic of the chosen Reynolds number (Re  = 150), is presented.
Time histories of lift (plot a) and drag (plot b) coefficients for 500 ..are reported in figure 8.The figure also includes a power spectral density (PSD) of the lift coefficient (plot c), with the Strouhal number St  =    ∞ related to the diameter of the three cylinders, and a representation of the trajectory in the time-delayed embedding space of the force coefficients (plot d).In the chaotic regime of the fluidic pinball, the main peak in the power spectral density of   is notably broad.However, a predominant one is found corresponding to St  = 0.148, associated with vortex shedding (with a shedding period  ℎ = 6.76 ..).The curve also shows a second peak associated with a lower Strouhal number (St  = 0.013) due to resonance in the flow field related to the finite size of the domain, as also observed in Deng et al. (2020).
The condition of flow symmetry ( C ≈ 0) is recovered in the chaotic regime, unlike the lower Reynolds number regimes which exhibit non-symmetric wakes (Deng et al. 2020).In addition,   standard deviation is (  ) = 0.112.As for the   , its free dynamics exhibits C = 3.46 and (  ) = 0.0658.Figure 9 shows the controlled case in ideal measurement conditions (i.e.no noise), with the control parameters automatically selected using the Bayesian optimization algorithm.The simulations presented from now on include an initial unforced phase before activating the control at time instant   = 50 ...The plots on the left show the time histories of the controlled aerodynamic coefficients and the exogenous input provided to the system 3 The control automatically selects an implementation law with the two rear cylinders in counter-rotation.The upper cylinder rotates clockwise ( 2 = − 3 < 0) with nearly constant rotational speed equal to the maximum values set in the optimization problem.This actuation mechanism corresponds to boat tailing, which has been previously observed and studied in various works on control of the wake of fluidic pinball (Raibaudo et al. 2020;Li et al. 2022).
This mechanism primarily aims to reduce the pressure drag of the fluidic pinball.Under unactuated conditions, as observed in plot  of figure 9, an extended recirculation region with low pressure and velocity forms behind the rear cylinders.The boat tailing redirects the shear layers of the top and bottom cylinder toward the pinball axis.This streamlining process of the wake increases the pressure behind the cylinders, delays the separation, and energizes the wake.Furthermore, the gap jet between the cylinders is substantially weakened.A detailed description of this mechanism can be found in Geropp & Odenthal (2000).A full wake stabilization is not achieved, i.e. the vortex shedding is not suppressed, due to the imposed constraints of   ≤ 1.This is in line with the study of Cornejo Maceda et al. (2021), where wake stabilization in boat tailing is achieved with rotation of the rear cylinder of 2.375 or larger.
A condition of global optimum in terms of minimizing the coefficient of drag would be achieved with an actuation that involves the rear cylinders in boat tailing and the front An examination of the behavior of the drag coefficient (  ) resulting from the implementation of the MPC control shows a decrease in its mean value, reducing it by 43.3% compared to the unforced case.Moreover, the standard deviation of the drag coefficient is reduced as well, by 85.8% during forced conditions compared to the unforced ones.
Figure 10 shows the out-of-plane vorticity in three different phases: before, during, and after the application of the control to the fluidic pinball wake.The MPC actuation is able to reduce the interaction between the upper and lower shear layers.In addition, the wake behind the fluidic pinball exhibits a more regular pattern in the post-transient phase.The wake meandering is strongly reduced, possibly due to the front stagnation point control with the front cylinder.The wake streamlining due to the boat tailing effect might also be contributing in this direction.This is also observed in the oscillating behavior of the   , which becomes more regular and with less intense peaks once the control is activated.Indeed, in the controlled phase, the power spectral density of the   presents a single very narrow peak concentrated at a Strouhal number of St  ≈ 0.14, resulting in a shedding period of 7.16 .., very similar to the predominant shedding frequency in the free response.
An analysis of the total power trend is required to assess if the MPC is able to achieve a net energy saving.The total power (  ) is characterized by the sum of two contributions: the first is directly related to the drag of the fluidic pinball (  ); the second one is associated with the power required for actuation (  ).Thus, the total power can be defined as follows: being   the torque acting on the   ℎ cylinder of the fluidic pinball.The time history of power during the free and forced phases is given in figure 11.It can be observed that under undisturbed conditions, the total power remains around an average value of 1.75, with oscillations related only to the drag of the fluidic pinball.When moving to the controlled solution, the total power undergoes a sharp decrease in the transient phase thanks to the boat tailing.The total power stabilizes at about 1.12, experiencing a reduction of 36.31% compared to the unforced phase.It is also observed that the average actuation power corresponds only to 12.13% of the total power in the controlled phase.The average value of the actuation power is mainly related to the rear cylinders in boat tailing, while its oscillations are due to the front cylinder aiming to control the shedding, in order to reduce the lift fluctuations.
It must be remarked though that this does not take into account possible inefficiencies that would arise in an implementation in a real environment.
A comparison between the proper orthogonal decomposition (POD) of the flow fields with and without control according to the MPC-based control algorithm is also shown in figure 12.It is specified that the POD has been performed on a dataset including both the transient and the stationary part of the control.On the left side of the figure, the spatial distributions of the 1st, 3rd, and 5th out-of-plane vorticity modes are shown for both free and controlled cases.The plots on the right side of the figure show the eigenvalues of the first 100 temporal modes.The eigenvalues are normalized with respect to their sum for the case without actuation.This allows a direct comparison of the mode energy in unforced and forced cases.
By observing the graphs related to the POD of the wake dataset under controlled conditions, it is noted that shedding inception is shifted upstream if compared to the unforced condition and is also more regular in the presence of boat tailing of the rear cylinders.Furthermore, the modes of the controlled case, compared to the free case, are characterized by lower energy, with the exception of the first 2 modes.This is ascribed to a partial wake stabilization, with reduced meandering.Mode 1 spotlights indeed more defined vortical structures, with crosswise width practically constant throughout the wake development.The fluidic boat tailing of the wake has indeed the effect of re-energizing the outer shear layers, redirecting them towards the pinball axis.The wake dynamics shift towards a less chaotic behavior, as also observed in the more regular oscillations of the   after activation of the control (figure 9).A less rich dynamics results in a more compact POD eigenspectrum for the controlled case, with higher energy in the first two modes.The reduced crosswise oscillations are observed also in the structure of higher-order modes.Modes 3 − 6 for the uncontrolled case embed the energy of higher-order harmonics and model the crosswise wake oscillations.In addition to the fluidic boat tailing due the the rotation of the rear cylinders, the stagnation point control enforced by the front cylinder has the effect of weakening such oscillations.This results in lower energy for the corresponding POD modes in the controlled case, and a more compact structure of their vortical features.

The effects of measurements noise in the sensors
The self-tuned control strategy with online LPR of the sensor signal is assessed in the presence of noise.The control results in terms of   and   and actuation are given in figure 13.We include as a reference the case of MPC without hyperparameter tuning, i.e. the control parameters are manually selected and all set to unity (except for the prediction/control window set to 3..).The outcomes when the parameter selection is automatically performed using BO is reported for different noise level.Without hyperparameter tuning, the performance of the control is quite poor for both   and   .In this case, the mean drag coefficient is 3.0598, i.e. significantly larger than the value of ≈ 1.97 achieved after hyperparameter tuning.Furthermore, MPC converges to an asymmetric controlled state, with a negative average lift coefficient.Quite surprisingly, this is achieved even if   2 =   3 = 1, i.e. forcing symmetry in the hyperparameters does not necessarily lead to symmetric actuation.After hyperparameter tuning, the penalty weight of actuation (and actuation variability) on the rear cylinder is decreased with respect to the front cylinder (see Table 1).This allows stronger actuation on the rear cylinders, thus fostering strong fluidic boat tailing, with consequent strong drag reduction, accompanied by a weak action of the front cylinder to improve wake stability.The cases analyzed above are reproduced also in figure 14 in terms of trajectories of   (),   (), and   ( − ) where  corresponds to a quarter of the shedding period  ℎ of the fluidic pinball wake.The plots offer a description of the attractor onto which the wake dynamics of the fluidic pinball evolves, including the transitional phase from free to forced conditions.When the parameter tuning is performed, during the control   experiences a remarkable reduction, whereas   evolves on a dynamics similar to the unforced case, although less chaotic.Nevertheless, it is worth noting that the control of   is not perfectly attained, as the model uncertainties hinder a complete description of the shedding dynamics of the fluidic pinball.Without hyperparameter tuning, in the first part of the transient the system transitions initially towards a state with a lower average   (around 2.5) and   oscillations of slightly lower intensity.The relatively high penalty on the actuation of the rear cylinders with respect to the front one, however, forces the system to avoid strong boat tailing and starts leveraging asymmetric actuation based on the front cylinder in the attempt to stabilize the wake through stagnation point control.The system finally converges to a limit cycle with a higher average   and with large oscillations of the   around a negative average.
In terms of the   , the control remains nearly unaffected by variations in measurement noise in the sensors.Regardless of the noise intensity, the posterior cylinders consistently maintain their boat tailing configuration, exhibiting maximum rotational velocity limited by the actuation constraints.On the other hand, the control of the lift coefficient demonstrates higher sensitivity to measurement noise.This is attributed to the increased susceptibility of the predictive model to disturbances in the initial conditions of online predictions.As the measurement noise intensifies, the frontal cylinder displays larger oscillations in the selected actuation law prescribed by the MPC-based control.
Nevertheless, LPR provides advantages to the algorithm, resulting in only marginal degradation of control performance.This increases the robustness of the application in realistic control scenarios characterized by high sensor noise.The panel in a) shows the results when the cost function parameters are hand-selected and equal to the unity, while panels b-d) show the result when Bayesian Optimization is used for hyperparameter tuning.Cases a) and b) relate to ideal measurement conditions while cases c) and d) have noise of 1% and 5%, respectively.LPR is included in the MPC control framework in all cases where measurement noise is present.

Symmetry in rear cylinder parameters
This part is intended to provide more details about the choice of symmetry in the control parameters.Specifically, it concerns the components of the weight matrices that penalize the input expenditures and input variability of the rear cylinders during control, already introduced in § 2.1.4.
The results of the control simulations with and without imposed symmetry conditions are presented in figure 15.A certain bias is observed in the control of   in the case without imposing symmetry conditions.This is related to an asymmetry in the rotation of the rear cylinders.However, under imposed symmetry, the control of the drag coefficient improves slightly, as the tangential speeds of rotation of the rear cylinders are both at the actuation limits, creating a stronger boat tailing.This lower effectiveness of the optimization in the case without imposed symmetry can be explained by the relatively low importance of the   fluctuations in the functional of Bayesian optimization in equation (2.10), as discussed in § 3.2.The implementation of a law that binds the parameters of the rear cylinders reduces the  total number of parameters to be tuned, leading to a less computationally expensive offline phase of the control.

Discussion and conclusion
In this work, a noise-robust MPC-based control algorithm that does not require user intervention neither for the plant modeling nor for the hyperparameters selection is proposed.The model used in the MPC optimization is selected from input-output data of the system under control using nonlinear system identification.The control parameters are automatically selected using a black-box optimization based on Bayesian methods.Additionally, the robustness of the algorithm is enhanced by the nonparametric smoothing technique LPR, which acts on the output data from the control sensors to address the presence of measurement noise.
The proposed control algorithm is successfully applied to the control of the wake of the fluidic pinball for drag reduction in a chaotic regime (specifically, for Re  = 150) using solely aerodynamic forces to guide the control strategy.The methodology achieves good success in controlling a non-linear, chaotic, and high-dimensional system.Being based on the MPC technique, the algorithm easily allows for the inclusion of non-linear constraints in the control and is very promising for applications where the control target is not a simple setpoint stabilization.
The proposed technique has shown to be highly robust to sensor measurement uncertainty, performing excellently even in realistic control scenarios characterized by a high level of noise.It must be remarked that, although rather realistic for force measurements, the Gaussian noise model might not be adequate if other sensing techniques are used.This may have an impact on plant model identification and LPR performance.However, the most attractive feature is that the method requires minimal user interaction, as the control parameters are automatically tuned by the Bayesian optimization algorithm according to the control target, taking into account also different fidelity levels for the plant predictions.The automatic procedure identified the most rewarding directions to optimize parameters with the aid only of a few numerical experiments.To control   and   of the fluidic pinball around a zero setpoint the MPC algorithm used a combination of two control mechanisms that have been previously considered in fluidic pinball control works: the boat tailing of the rear cylinders for drag reduction and the stagnation point control of the front cylinder for shedding control.This corresponded to a stronger penalty on the   and a low penalty on the actuation cost of the rear cylinders.Without the framework proposed in this work, this process would have required parametric studies or sub-optimal analysis, with the inherent difficulty of choosing relative weights between heterogeneous quantities.It is worth highlighting that the control algorithm is currently being tested on a control case with relatively simple dynamics (Deng et al. 2020(Deng et al. , 2022)).This streamlines the process of coordinate selection and nonlinear system identification using the SINDYc technique.On the one hand, the aerodynamic force coefficients in separated wake flows exhibit features that render them suitable for plant modeling.On the other hand, the dynamics of the flow, although chaotic, evolve clearly on a low-dimensional attractor, thus simplifying the library selection for the application of SINDYc.For more complex flows (higher Reynolds number or with less clear features for straightforward coordinate identification) we envision that the challenge of MPC application will reside mostly in the coordinate selection and in the plant identification.Fast-paced advancements in grey-and black-box modeling are paving the way to interesting research pathways in these directions.The approach proposed here, however, will still be applicable and will benefit from such advancements.The offline optimization of the hyperparameters of the cost function can be applied independently of the method for coordinate selection or plant identification.Furthermore, the outcome will automatically adapt to the fidelity of the model plant when compared to the uncertainty of the measurements used to feed it.
Funding.The authors acknowledge the support from the research project PREDATOR-CM-UC3M.This project has been funded by the call "Estímulo a la Investigación de Jóvenes Doctores/as" within the frame of the Convenio Plurianual CM-UC3M and the V PRICIT (V Regional Plan for Scientific Research and Technological Innovation).
Figure 16: MPC-based control algorithm schematic.Steps for implementing the control algorithm are shown, including dataset generation, creation of the predictive model using SINDYc, and parameter tuning.The main block displays the closed-loop MPC scheme, including also LPR for the reduction of the effects of sensor measurement noise.

.
The authors thank Prof. B. Noack, from Harbin Institute of Technology, and Prof. M. Morzynski, from Poznan University of Technology, for providing the code for the fluidic pinball simulations, and Dr. G. Cornejo Maceda for the post-processing code for force estimation.

Figure 1 :
Figure 1: MPC-based control algorithm schematic: dataset generation, creation of the predictive model and parameter tuning.The main block displays the closed-loop MPC scheme, including also LPR to mitigate the effects of sensor measurement noise.

Figure 17 :
Figure 17: Domain of the incompressible two-dimensional DNS of the flow around fluidic pinball.Front, top and bottom cylinders are respectively labeled as 1, 2, and 3.The rotational velocities of the cylinders are  1 ,  2 , and  3 , respectively.The arrows indicate positive (counterclockwise) rotations.The background shows the 8k nodes grid used for the DNS.The contour colors indicate the out-of-plane vorticity.

Figure 3 :
Figure 3: Domain of the incompressible two-dimensional DNS of the flow around fluidic pinball.Front, top and bottom cylinders are labeled as 1, 2, and 3, respectively.The rotational velocities of the cylinders are  1 ,  2 , and  3 .The arrows indicate positive (counterclockwise) rotations.The background shows the 8633 nodes grid used for the DNS.The contour colors indicate the out-of-plane vorticity.

Figure 4 :
Figure4: Prediction performance of the force model obtained using SINDYc.The plots display the average value and the confidence region for  of the probability distribution of the normalized   and   prediction errors with respect to their standard deviation under unforced conditions.Panels on the left show the results in absence of measurement noise of the initial condition while those on the right with 1% of measurement noise.In the presence of noise, a comparison is shown with the results of predictions obtained using LPR estimation as initial condition.

Figure 5 :Figure 6 :
Figure5: LPR applied to a   time series within a time period of 60...This case corresponds to forced fluidic pinball dynamics.Top panels display   estimation in the presence of increasing noise levels (1%, 3%, 5%).Bottom panels show LPR estimation of   .Each plot includes the ideal and noisy time series, also including the RMSE of the LPR estimation.

Figure 7 :
Figure 7: Representation of the cost functional   in equation (2.10), the component due to   control ( 2  ), and the component due to   control ( 1  ).Contour plots are presented as functions of  1 and  2 (top graphs) and   1 and   2 =   3 =   2,3 (bottom graphs).The control case considers a 1% measurement noise and symmetric weighting coefficients for the rear cylinders of the fluidic pinball.In the graphs, the remaining components of the hyperparameter vector  are fixed at their optimal values.

Figure 8 :
Figure 8: Characteristics of the flow field around the fluidic pinball in undisturbed conditions.Panels a-b) show the global lift and drag coefficients of the fluidic pinball.Plot c) shows a PSD of the lift coefficient while d) provides a representation of the trajectory in the time-delayed embedding space of the force coefficients   (),   () and   ( − ), being  a quarter of the shedding period of the fluidic pinball wake.The peak in the PSD is at a Strouhal number St  = 0.148.

Figure 9 :
Figure 9: Results of the control application in absence of measurement noise in the sensors.Graphics on the left show the time histories of the   ,   and input vector  during an initial unforced phase and then during an active control phase.Graphics on the right show the mean value of the streamwise component of the velocity () in the time frames denoted A, B and C, respectively.Black arrows indicate the direction of rotation of the cylinders.

Figure 10 :
Figure 10: Contour plot of the out-of-plane vorticity component of the flow around the fluidic pinball.Several instants corresponding to the phases of control onset, transient, and post-transient phase are presented. ℎ refers to the characteristic shedding period of the fluidic pinball in unforced conditions.

Figure 11 :
Figure11: Time histories of the power associated with drag, actuation, and total power during the free and forced stages.Simulation of the fluidic pinball wake control in the absence of measurement noise.

Figure 16 :
Figure 16: Results of a POD applied to the flow field of the wake past fluidic pinball.Panels on the left depict the representation of out-of-plane vorticity for the spatial modes 1 (plots a), 3 (plots b) and 5 (plots c) under both free and forced conditions in accordance with the MPC-based control architecture.Every plot also presents velocity vectors corresponding to each spatial mode.The graph on the right shows the squared singular values (  ) of the POD, normalized with respect to the sum of the squared singular values for the free case ( ,    ).POD performed on a dataset consisting of  = 1200 snapshots of the fluidic pinball wake, including the transient in the forced case.

Figure 12 :
Figure 12: Results of a POD applied to the flow field of the wake past fluidic pinball.Panels on the left, from top to bottom depict the representation of out-of-plane vorticity for the spatial modes 1, 3 and 5 under both free and forced conditions in accordance with the MPC-based control architecture.Every plot also presents velocity vectors corresponding to each spatial mode.The graph on the right shows the squared singular values (  ) of the POD, normalized with respect to the sum of the squared singular values for the free case ( ,    ).POD performed on a dataset consisting of  = 1200 snapshots of the fluidic pinball wake, including the transient in the forced case.

Figure 1 :Figure 13 :
Figure 1: Time series of   ,   and exogenous input   at free and forced stage according to the MPC based control framework applied to the fluidic Pinball.The panels in a) show the results when the cost function parameters are hand-selected and equal to the unity.The panels in b) show the result when Bayesian Optimization is used for hyperparameter tuning.LPR is included in the MPC control framework in all cases where measurement noise is present.

Figure 14 :
Figure14: Trajectories in the time-delayed embedding space of the force coefficients   (),   () and   ( − ), being  a quarter of the shedding period of the fluidic pinball wake.The panel in a) shows the results when the cost function parameters are hand-selected and equal to the unity, while panels b-d) show the result when Bayesian Optimization is used for hyperparameter tuning.Cases a) and b) relate to ideal measurement conditions while cases c) and d) have noise of 1% and 5%, respectively.LPR is included in the MPC control framework in all cases where measurement noise is present.

Figure 15 :
Figure15: Results of the application of MPC-based control in terms of   ,   , and input vector , without (left panels) and with (right panels) the imposition of symmetry in the parameters related to the rear cylinders.In both cases, all the parameters are tuned using BO algorithm.Results of control in the absence of measurement noise in the sensors.
The control window   is shown in orange.Dashed lines indicate future state and actuation predictions.Blue circles represent a discrete sampling of the system state prediction.The continuous formulation allows nonmandatory discrete sampling and step-like actuation can be relaxed.
() ′ ,  ∈ B ⊂ R   , where B is the set of allowable inputs.Denoting the time derivative of the state vector with , the system p Figure 2: Graphical representation of MPC strategy for stabilizing around a set-point (horizontal dashed line).Past measurements (light-blue-shaded region) depict system state (blue lines with squares) and actuation (green line).

Table 1 :
.10) and selecting only Results of the control tuning through Bayesian optimization.A comparison is