Pierre Clément Blaud

Philippe Chevrel

Fabien Claveau

Pierrick Haurant

Anthony Mouraud

From multi-physics models to neural network for predictive control synthesis

Keywords: Feedforward neural network, hyperparameters optimisation, economic model predictive control, modelica, model identification

Introduction

Systems and control technology are everywhere, from homes to airplanes, and they face great challenges in areas, such as water networks, increase share of renewable energy, intelligent transport, process engineering or medical science. Some challenges require interdisciplinary methods ranging from advanced control techniques, modelling interconnected complex systems and machine learning [START_REF] Lamnabhi-Lagarrigue | Systems & Control for the future of humanity, research agenda: Current and future roles, impact and grand challenges[END_REF]. Among advanced control techniques, Model Predictive Control (MPC) is commonly studied in academia [START_REF] Mayne | Model predictive control: recent developments and future promise[END_REF] and used in industrial applications [START_REF] Darby | MPC: current practice and challenges[END_REF], because it enables considering constraints on the state variables of the process. Since initial implementation in the petroleum industry [START_REF] Richalet | Model predictive heuristic control: application to industrial processes[END_REF], MPC has broadened its field of industrial applications to include traffic control [START_REF] Ye | A Survey of Model Predictive Control Methods for Traffic Signal Control[END_REF], energy management in residential house [START_REF] Lefort | Hierarchical control method applied to energy management of a residential house[END_REF], heating ventilation and air-conditioning control [START_REF] Afram | Theory and applications of HVAC control systems -A review of model predictive control (MPC)[END_REF], aerospace [START_REF] Weiss | Model predictive control for spacecraft rendezvous and docking: strategies for handling constraints and case studies[END_REF] or the automotive industry [START_REF] Hrovat | The development of model predictive control in automotive industry: a survey[END_REF]. Classical MPC involves a state space model of the dynamical system, states constraints, inputs constraints and a quadratic cost function to minimize. The optimization is performed online and the first optimal control sample is the one applied to the actuators, while a new optimization is relaunched to prepare the next step. There are some extensions or variants of MPC, such as Economic Model Predictive control (EMPC), where an economic cost is substituted with the quadratic one in order to control the dynamical system for economic behavior [START_REF] Rawlings | Fundamentals of economic model predictive control[END_REF][START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF].

Finding a suitable dynamical model allowing predictive control synthesis poses a challenging task for control engineering. Models that finely represent the process are generally complicated, which makes it difficult to use them for MPC, and model simplification is often required. To acquire a suitable model, one method is to linearize the physical model around an operative point, which allows controller tuning and leads to linear MPC [START_REF] Muske | Model predictive control with linear models[END_REF] or linear EMPC [START_REF] Jørgensen | Economic MPC for a linear stochastic system of energy units[END_REF]; however, the model poorly represents the non-linear dynamical system. To increase model accuracy, a linear time-varying model from the non-linear model can be adopted [START_REF] Murillo | Iterated non-linear model predictive control based on tubes and contractive constraints[END_REF] or more directly the non-linear model [START_REF] Grüne | Nonlinear model predictive control[END_REF]. Furthermore, a hybrid MPC is possible to handle systems with continuous states, logic rules, discrete states or if-then-else conditions [START_REF] Lazar | Model predictive control of hybrid systems : stability and robustness[END_REF]. The engineering time required for the modelling process of modelling is labour intensive and may be out of proportion based on industrial expectations [START_REF] Maddalena | Data-driven methods for building control -A review and promising future directions[END_REF].

One method of reducing engineering time is to identify the model from inputoutput signals of the process to achieve a black-box model. Black-box models used for predictive controllers can be canonically parametrized model or a neural network. An AutoRegressive Exogenous (ARX) model is presented in [START_REF] Huusom | Tuning of methods for offset free MPC based on ARX model representations[END_REF], while a recurrent neural network is used by an MPC in [START_REF] Wong | Recurrent neural network-based model predictive control for continuous pharmaceutical manufacturing[END_REF] and also in [START_REF] Wu | Economic Machine-Learning-Based Predictive Control of Nonlinear Systems[END_REF] with an EMPC. In [START_REF] Terzi | Learning-based predictive control of the cooling system of a large business centre[END_REF], the authors compared parametric models and neural networks models to identify a cooling system for a business center and found that the neural network model is more accurate compared to a parametric one. In [START_REF] Wong | Recurrent neural network-based model predictive control for continuous pharmaceutical manufacturing[END_REF], the data were generated from the non-linear state space model, by simulation. In [START_REF] Terzi | Learning-based predictive control of the cooling system of a large business centre[END_REF], the data originate from experimental measures and yield favorable results; however, there are many cases where measured data are lacking or expensive [START_REF] Schrangl | Iterative Model Identification of Nonlinear Systems of Unknown Structure: Systematic Data-Based Modeling Utilizing Design of Experiments[END_REF].

The method considered in this work is intended for dynamical systems to be controlled, for which measured data are lacking, and therefore data will be generated via simulation using multi-physical modeling. Model identification is then performed to acquire a model which requires less calculation time for its simulation while maintaining accuracy. The following method includes several advantages: First, when using a simulation model, a dedicated researcher or engineer performs physical modelling using specific simulation tools, which makes the simulation model accurate regarding the state of the art. Second, this method reduces engineering time for the control researchers in finding a blackbox model rather than a classical physical model for a predictive controller. Third, this method helps finding a black-box model with greater accuracy than classical physical models. Fourth, in some real systems, there are insufficient measurement data to identify an accurate black-box model, this method helps in generating data. Moreover, many simulation tools represent dynamical system in detail, such as the Modelica language [START_REF] Mattsson | Physical system modeling with Modelica[END_REF], which allows modeling multiphysical networked systems such as mechanical, electrical, thermal, hydraulic, pneumatic and fluid. In addition, commercial or open-source libraries model many systems, such as electrical [START_REF] Franke | Flexible modeling of electrical power systems-the Modelica PowerSystems library[END_REF], buildings [START_REF] Wetter | Modelica buildings library[END_REF], power plant system [START_REF] Hefni | Dynamic modelling of a combined cycle power plant with ThermoSysPro[END_REF], district heating system [START_REF] Giraud | Presentation, validation and application of the District Heating modelica library[END_REF] or greenhouses [START_REF] Altes-Buch | Greenhouses: a modelica library for the simulation of greenhouse climate and energy systems[END_REF].

Despite the popularity of i) multi-physics modeling for simulation on the one hand, ii) neural network modeling and iii) predictive control method, the current literature seemed to the authors not to make the most of the recent advances on these subjects, in particular by not proposing an integrated methodology for the global synthesis of control laws. The main contribution of this work comes in the form of a methodology systematically addressing the issue of predictive controller synthesis relying on a dynamic model based on a neural network for its implementation, this model being identified from data from multi-physical simulations (in the absence of big measurement data). The identification of the process proceeds from hyperparameters optimized with a meta-heuristic algorithm to obtain the ad hoc neural network. The predictive controller thus deduced from the application of the proposed methodology was successfully applied to a conventional control benchmark. In summary, the methodology proposed in this article links:

• The use of multi-physics modeling by taking advantage of its ability to combine physical subsystems of different types as well as the capitalization and reuse of such models thanks to shared libraries;

• The possibilities of identifying neural networks based non-linear dynamic models, with an automated choice of hyperparameters, the latter being heuristically optimized according to explicit criteria;

• The possibilities offered in control theory by non-linear model-based predictive control;

The paper is organized as follows: in Section 2, main elements of the methodology are depicted; in Section 3, the methodology is detailed; Section 4 presents the application of the methodology to a classical control benchmark; the results are depicted in Section 5; Section 6 concludes the work.

Methodology highlights

This study's methodology is shown in Fig. 1 and divided into five stages. This section provides a tutorial to aid understanding and implementation. The first step regards the fine-scale modelling using multi-physical modelling tools, which can be performed by modelling experts to achieve a model with the closest resemblance. In addition, the model's parameters can be scaled using known completed plant parameters and using calibration from plant data measures if some are missing. Information about complex networked physical systems modeling is provided in Section 3.1.

Data generation from simulation results

System

The second step involves data generated from the simulation model designed in the first step. The model is extensively simulated from selected scenarios to best cover the conditions of the system's use to be controlled, and the data produced must be both qualitative and quantitative to obtain a representative neural network (by optimization). These aspects are discussed in Section 3.2.

The third step concerns identifying the neural network approximating the System x(t) u(t) dynamic behavior of the process (see Section 3.3). The accuracy of the targeted predictive control depends on the quality of this model, which will be used as an internal model. The description of this model and methods analyzing its stability will be specified, and an automated tuning method of an appropriate set of hyperparameters will be introduced to aid for building and optimizing the neural network.

The fourth step regards with establishing the predictive control. In this work, the predictive control considered is economic, although the methodology could be applied to a classical MPC with a quadratic cost function. Section 3.4 introduces the EMPC (cost and constraint function) and specifies using the neural model from step three. In addition, some computer tools for implementing the controller and resolving the optimization problem will be detailed.

All steps are detailed in the following sections.

Methodology

Modelling complex systems

The modelling of complex networked physical systems has been a focus for years through formalism as bond graph [START_REF] Gawthrop | Bond-graph modeling[END_REF] or Energetic Macroscopic Representation (EMR) [START_REF] Agbli | Multiphysics simulation of a PEM electrolyser: Energetic Macroscopic Representation approach[END_REF] that enable graphical representation. Both are based on energy flows, effort and power relation, and physical relations use time-dependent Ordinary Differential Equation (ODE) supporting 0-D modelling. In addition, a computing language such as Modelica allows modeling and simulation based on a uniform general language for model design. A graphical representation of each sub-system and their interactions may be used, either from bond graph [START_REF] Vasilyev | Component-based modelling of PEM fuel cells with bond graphs[END_REF] or Modelica's graphical representations. Each block and subsystem are assembled (graphically or not) to build the whole system [START_REF] Mattsson | Physical system modeling with Modelica[END_REF]. A Modelica code concatenating all the physical equations is produced, which enables numerically simulation using a numerical solver such as a Differential Algebraic System Solver (DASSL) [START_REF] Liu | Comparison of 4 numerical solvers for stiff and hybrid systems simulation[END_REF].

Data generation

Generating simulation data is mandatory for system identification. Before generating data, the control actuators must first be selected and will form the inputs of the identified dynamical system; second, the measures must be selected and will form the states of the identified dynamical system. As a result, we will use an input-output description of the system (see Fig. 2). Both actuators and measures must be chosen wisely according to the dynamical system capabilities and with physical implementation possibilities. After inputs and outputs are chosen, the signals sequences are applied to the inputs while the outputs are measured. The characteristics of input sequences of the simulated dynamical system require careful selection to generate output data that reveal its dynamic nature, which is called persistent excitation [START_REF] Schrangl | Iterative Model Identification of Nonlinear Systems of Unknown Structure: Systematic Data-Based Modeling Utilizing Design of Experiments[END_REF]. Various inputs could be used to excite the system, such as constant input signals (on-off), ramps, steps, sinusoidal, random sequences or Pseudo-Random Binary Sequences (PRBS) [START_REF] Schrangl | Iterative Model Identification of Nonlinear Systems of Unknown Structure: Systematic Data-Based Modeling Utilizing Design of Experiments[END_REF][START_REF] Landau | System identification: the bases[END_REF]. These sequences lack the same persistent excitation. In [START_REF] Schrangl | Iterative Model Identification of Nonlinear Systems of Unknown Structure: Systematic Data-Based Modeling Utilizing Design of Experiments[END_REF], the authors found that random sequences adequately excite the system; in [START_REF] Landau | System identification: the bases[END_REF], the authors indicate that PRBS can excite the dynamical system with distinct sinusoidal components. In addition, a sequence that is a mix of the random sequence and PRBS is called a Multilevel Pseudo-Random Signal (MPRS).

Identification 3.3.1. Dynamical system

In this work, we considered a discrete time invariant system of the following form:

xk+1 = f (x k , ūk ) (1)
where x ∈ Xnx is the system state, ū ∈ Ūnū is the command system input and k ∈ I ≥0 is the discrete sampling time occurrences. Moreover, all states are available for the controller.

Feedforward neural network

The first mathematical model of a neuron was proposed by W. McCulloch and W. Pitts in 1943 and is called the McCulloch-Pitts model [START_REF] Yang | Nature-inspired metaheuristic algorithms[END_REF]. The mathematical formulation of single artificial neuron with multi-inputs is equal to [START_REF] Yang | Nature-inspired metaheuristic algorithms[END_REF]:

a = σ(W p + b) ( 2
)
where a is the output of the neuron, σ is the activation function, W is the weight matrix, p is the input of the neuron and b is the bias. One neuron is clearly insufficient to extract interesting characteristics from the data, which are usually stacked and leads to a Feedforward Neural Network (FNN) [START_REF] Kim | Standard representation and unified stability analysis for dynamic artificial neural network models[END_REF]. A schematic representation of FNN is shown in Fig. 3. FNNs represent the simplest neural network architectures and were chosen for implementation in this work. The mathematical formulation of a FNN with i hidden layers with the purpose to identify the dynamical system from Eq. ( 1) is equal to:

x k+1 = f h1 (v k ) • f h2 (a 1 ) • • • • • f hi-1 (a i-2 ) • f hi (a i-1 ) (3)
where x k+1 is the predicted output of the dynamical system by the neural network, and f hi is the hidden layer, such as

a i = f hi (p i ) = σ i (W i p i + b i
) with a i the output and p i the input. In addition, σ i is the activation function, W i is the weighting matrix and b i is the bias. Then, v k is the input of the FNN, which is the concatenation of the dynamical system inputs ūk and states xk at current and past points in time, v k = cat(x k , ūk , xk-1 , ūk-1 , ..., xk-l , ūk-l ), with l the back horizon as number of lag.

x k+1 u k x k u k-l x k-l ⋮ Figure 3
: Schematic view of a feedforward neural network with four hidden layers, inputs layer neurons are black, hidden layers neurons are grey, and output layer is white.

Fidelity measure

The loss function is employed to measure the quality of the FNN to reproduce a signal. This provides a quantitative scoring that shows the degree of similarity between from the data and the FNN output [START_REF] Wang | Mean squared error: Love it or leave it? A new look at Signal Fidelity Measures[END_REF]. Several loss functions are available, such as the Mean Square Error (MSE):

f loss = 1 N b N b i=1 [x i -xi ] 2 (4)
where N b is the number of samples, x i is the neural network output and xi is the target. Other loss functions include Mean Absolute Error (MAE), Root Mean Square Error (RMSE) and Mean Absolute Percentage Error (MAPE) [START_REF] Li | On comparing three artificial neural networks for wind speed forecasting[END_REF].

Hyperparameters tuning

The hyperparameters defining the neural network influence its performance and include the number of hidden layers, number of neurons per hidden layers, type of activation function, batch size, epochs, the optimizer used and its own parameters [START_REF] Yang | On hyperparameter optimization of machine learning algorithms: Theory and practice[END_REF]. For widelyd use accelerated gradient descent optimizers, parameters include the learning rate and momentums [START_REF] Zhang | Yellowfin and the art of momentum tuning[END_REF]. It is difficult to tune hyperparameters since they offer numerous choices, and we propose to automate these choices by selecting them via optimization, therefore at a higher hierarchical level, to minimize the final loss function. Performing optimization is challenging since it is non-differentiable and involves mixed and constrained variables. The optimization problem is formalized in the form:

min χ f loss (χ) (5a) subject to f nn (σ, δ, θ, µ, ν, o, τ, υ, φ) (5b) σ ∈ Σ (5c) δ ∈ ∆ (5d) θ ∈ Θ (5e) µ ∈ M (5f) ν ∈ N (5g) o ∈ O (5h) τ ∈ T (5i) υ ∈ Υ (5j) φ ∈ Φ (5k)
where f loss is the loss function for a given data set. f nn is the neural network architecture with dedicated inputs and outputs and Eq. (5b) is the related constraint. σ is the activation function and Σ is the constraint, see Eq. (5c). δ is the number of neurons per hidden layer and ∆ is the related constraint, see Eq. (5d). θ is the number of hidden layer and Θ is the hidden layer constraint, see Eq. (5e). µ is the epoch number and M is the constraint, see Eq. (5f). ν is the batch size and N is the constraint, see Eq. (5g). o is the optimizer and O is the constraint, see Eq. (5h). τ is the optimizer learning rate and T is the related constraint, see Eq. (5i). υ is the optimizer momentum exponential decay and Υ is the constraint, see Eq. (5j). φ is the optimizer momentum estimate and Φ is the related constraint, see Eq. (5k). Finally χ is the vector of optimization variables.

Implementation on computer machines

FNNs are implemented within machine learning tools such as PyTorch [START_REF] Paszke | Pytorch: an imperative style, high-performance deep learning library[END_REF] or Flux [START_REF] Innes | Flux: elegant machine learning with Julia[END_REF], which help reduce implementation complexity and allow training neural networks with Graphical Processing Unit (GPU) to accelerate training. In addition, a meta-heuristic algorithm is used to solve the optimization problem Eq. ( 5) due to the problem's of non-differentiability. The era of meta-heuristics features abundant algorithms for finding the best set of hyper-parameters, such as genetic algorithms or particle swarm optimization [START_REF] Abdel-Basset | Metaheuristic Algorithms: A Comprehensive Review[END_REF]. In addition, metaheuristic algorithms are available in language tools such as DEAP [START_REF] Fortin | DEAP: Evolutionary algorithms made easy[END_REF] or Black-BoxOptim.jl [START_REF] Feldt | BlackBoxOptim.jl[END_REF]. The implementation proposed in this paper is discussed in Section 4.2.

Economic Model Predictive Control 3.4.1. Description

There are two types of EMPC: The first involves a terminal cost and/or terminal states constraints [START_REF] Ellis | A tutorial review of economic model predictive control methods[END_REF], while the second has no terminal cost and constraints [START_REF] Grüne | Economic receding horizon control without terminal constraints[END_REF] and is employed in this work. The cost function used is moreover not quadratic, but associated with a so-called "economic" cost (Economic MPC). Traditionally, approaches seeking to guarantee the stability of MPC controls, introduce a constraint or a terminal cost. Traditionally, approaches seeking to ensure the stability of MPC controls introduce a constraint or a terminal cost. Here, a contractive constraint was added on the terminal state to ensure the MPC to be stabilized [START_REF] De Oliveira Kothare | Contractive model predictive control for constrained nonlinear systems[END_REF]. The formulation is described below, but please note that neural network input v k is reworded with x k , u k for notation simplicity :

min x,u N -1 k=0 l e (x k , u k ) (6a) subject to x k+1 = f * nn (x k , u k ) ( 6b
)
x 0 = x(t) (6c) x ∈ X (6d) u ∈ U (6e) x 0 P x0 ≤ x N P xN (6f)
where Eq. ( 6a) is the economic cost to minimize and N is the number of samples defining the length of the considered time horizon. Eq. (6b) is the time evolution constraint of the state, which is based on the best neural network architecture found in Section 3.3.4. Let's note that optimization problem ( 6) is non-linear, in particular, Eq. ( 6b) is non-linear due to FNN's activation functions. Moreover, we opted for a multiple shots approach. So, it is clear that the EMPC optimization stage is nonlinear, hence generically non-convex. It possibly leads to a only local optimum [START_REF] Grüne | Numerical Optimal Control of Nonlinear Systems[END_REF]. Eq. (6c) is the plant state measure, Eq. (6d) is the state constraints and Eq. (6e) is the input constraints. Eq. (6f) is the contractive constraint, P the contractive weighting matrix, and x is the state deviation from the optimal steady state. The EMPC is computed at each iteration, and the first sample of the optimal input computed is applied to the plant's actuators until the next sampling period. In addition, the EMPC defines an implicit feedback law:

µ N (x) = u * 0 (7)
Where µ N (x) is the implicit feedback related to x and u * 0 is the optimal input computed.

Computing machine implementation

The EMPC is translated and implemented in a modeling language specific to the field of optimization, and the optimization problem must be solved by an ad hoc solver. Several tools are available, such as JuMP [START_REF] Dumming | JuMP: a modeling language for mathematical optimization[END_REF] or Casadi [START_REF] Andersson | CasADi: a software framework for nonlinear optimization and optimal control[END_REF]. The neural network used to define the state evolution constraint is a non-linear model, and thus a non-linear programming solver is needed, such as the Ipopt solver, which implements the interior point method [START_REF] Wächter | On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming[END_REF]. Implementing such an on-line implicit feedback law µ N (x) with high sampling can be computationally difficult [START_REF] Abughalieh | A Survey of Parallel Implementations for Model Predictive Control[END_REF], especially when the model is non-linear [START_REF] Diehl | Efficient Numerical Methods for Nonlinear MPC and Moving Horizon Estimation[END_REF]. The optimization

Lower tank q a q b V a V b Tank 3 Tank 4
Tank 1 Tank 2 computation cannot exceed a sampling period. For some problems, a parallel implementation within computing machines composed of multi-core Central Processing Unit (CPU), GPU or Field-Programmable Gate Array (FPGA) [START_REF] Abughalieh | A Survey of Parallel Implementations for Model Predictive Control[END_REF] is necessary.

Example

Four tank benchmark

To illustrate the relevance of the methodology, we apply it to a benchmark widely used for ordering: the level control of a Quadruple Tank Process (QTP), also known as the four-tank benchmark [START_REF] Johansson | The quadruple-tank process: a multivariable laboratory process with an adjustable zero[END_REF], is widely used by academics [START_REF] Johansson | Teaching multivariable control using the quadruple-tank process[END_REF] and recognized for its capacity to illustrate industrial issues [START_REF] Shneiderman | Properties and control of the quadrupletank process with multivariable dead-times[END_REF]. Contributions include the use of a decentralized PID controller [START_REF] Hušek | Decentralized PI Controller Design Based on Phase Margin Specifications[END_REF] or a control by sliding mode in [START_REF] Prusty | Sliding mode control of coupled tank systems using conditional integrators[END_REF]. The implementation of a predictive control applied to the QTP benchmark is also widely discussed, and several distributed MPC methods are compared in [START_REF] Alvarado | A comparative analysis of distributed MPC techniques applied to the HD-MPC four-tank benchmark[END_REF]. The QTP benchmark has four tanks, two three-way valves and two pumps. The aim is to maintain a certain water level in the reservoirs. Pump q a fills tanks 1 and 4, then pump q b fills tanks 2 and 3. In addition, water leaks through tanks orifices, as tank 3 leaks to tank 1 and tank 4 to tank 2, then tanks 1 and 2 to lower tank, Fig. 4. As a result, the pump flows must be greater than zero. The main parameters of the benchmark are shown in Tab. 1.

Application of the proposed methodology 4.2.1. Physical modelling

The QTP is firstly modelled using Modelica to support the simulation tests, which is fine modeling that considers water property, pumps inertia, pipe roughness and valve pressure loss. Each sub-model is graphically assembled graphically, see Fig. 5. q min or q max m 3 /h

Data generation

The actuators of the dynamical system are the two pumps, and the measures are the four tanks water levels. Each water level in the reservoir can be viewed as a state variable: xi = h i . The inputs are the pump water flows ūi = q i . Two random sequences were applied to the control inputs: The first is a PRBS signal, whose properties are provided in Tab. 2; the second is the piecewise constant with constant random values ranging from minimum to maximum water flow, see Tab. 3. The simulations were performed to generate a total of 10,368,000 samples for each neural network input with a 5s sample period. Furthermore, the data are filtered to only use data within water level constraints, which leads to 5,686,738 data samples available to feed the neural network.

Dynamical system identification

The neural network inputs are water flows from pumps, and current water levels with back horizon l = 0, put of the neural network is the water levels predicted for the next sample time, x k+1 = hcat(h 1 k+1 , h 2 k+1 , h 3 k+1 , h 4 k+1 ). The data generated are separated into two arrays for the training stage, where the first array contains the neural network inputs while the second has the outputs. The available data [START_REF] Altes-Buch | Greenhouses: a modelica library for the simulation of greenhouse climate and energy systems[END_REF]686,738) are separated and randomized into 80% for training and 20% for testing, and during each of the 20 epochs, 10% of the training data are extracted for validation while the remainder feeds the neural network. Furthermore, the chosen loss function chosen is the MSE, and the neural networks were implemented with Flux.jl [START_REF] Innes | Flux: elegant machine learning with Julia[END_REF]. The meta-heuristic algorithm chosen to optimize the hyperparameters is a differential evolution algorithm [START_REF] Das | Recent advances in differential evolution-an updated survey[END_REF] with adaptive setting [START_REF] Al-Dabbagh | Algorithmic design issues in adaptive differential evolution schemes: Review and taxonomy[END_REF] q min to q max m 3 /h and radius limited from Julia package BlackBoxOptim.jl [START_REF] Feldt | BlackBoxOptim.jl[END_REF].Hyperparameters range values are shown in Tab. 4. Otherwise, the neural network training were performed on Nvidia P6000 GPUs thanks to JuliaGPU [START_REF] Besard | Effective extensible programming: unleashing Julia on GPUs[END_REF].

v k = hcat(q a k , q b k , h 1 k , h 2 k , h 3 k , h 4 k ).

EMPC tuning and implementation

EMPC is optimized every 5s and over a time horizon equal to 25s (5 periods). The economic cost considered is a function of the squared weighted sum of the pump flow rates and the inverse of the water volumes in reservoirs one and two [START_REF]A robust economic MPC for changing economic criterion[END_REF]. This reflects the objective of the search for two conflicting objectives: the maximization of the volume of water and the minimization of the energy consumption of the pumps (linked to the square of the flow rates):

l e (x k , u k ) = u 2 1 k + κ 1 u 2 2 k + κ 2 V min S c (x 1 k + x 2 k ) ( 8
)
Where κ 1 and κ 2 are time varying costs. V min = 0.012m 3 is the minimum volume of water in a tank, and S c is the tank cross-section. κ i = (κ 1 ; κ 2 ) values are equal to κ i = (1; 500) from 1 to 1,500s, κ i = (1; 1000) from 1,500 to 3,000s, κ i = (5; 100) from 3,000 to 4,500s and κ i = (0.5; 500) from 4,500 to 6,000s, Tab. 5. Also, the P matrix selected is the identity matrix and the optimal steady state are visible in Tab. 5. Furthermore, constraints enable considering the water tanks' max level (state constraints) and admissible water flows (input constraints), Tab. 1. The predictive control is implemented within JuMP.jl with multiple shooting numerical method and the solver used is Ipopt within a workstation equipped with an Intel Xeon Gold 5122 CPU.

Simulation test bench

The QTP fine-scale modelling from Modelica was exported using a Functional Mock-up Interface (FMI). A simulation software (Simulink) simulated the Functional Mock-up Unit (FMU) and saved states and inputs in text files throughout simulation. Furthermore, every 5s, the simulation was paused to execute the predictive control. Finally, the optimal control inputs computed were sent to the FMU and the simulation run again until the next iteration. A schematic representation is provided in Fig. 6.

Results

System identification

Fig. 7.a depicts the boxplot hyperparameters optimization results with validation loss versus trained networks, where 618 neural networks were trained, 487 algorithm steps were performed and vanishing gradients occurred with 18 trained networks. The lowest validation loss is equal to 8.11×10 -7 and the highest is equal to 3.79 × 10 -1 . The validation loss median is equal to 5.17 × 10 -6 , and as a result, 50% of the trained networks have validation loss higher than 5.17 × 10 -6 . Hyperparameters optimization allows selecting the best set of hyperparameters for a network with accurate identification. Fig. 7.b depicts the validation loss with clusters of 10 trained, which shows three groups: The first is from the first trained network until 200 trained networks. The boxplot validation loss slightly decreased and reached a lower validation loss at 200 trained networks. The second group, from 210 trained networks to 470 trained networks, validation loss increased and then decreased and reached lower validation loss at 470 trained networks. The third group, from 480 trained networks until the last trained networks, validation loss increased and then decreased and reached lower validation loss at 600 trained networks. This behavior shows that the meta-heuristic algorithm can select hyperparameters to reduce validation loss, but it needed a couple iterations. In addition, it reached some local minimal and tried to escape to find a lower validation loss. Furthermore, due to metaheuristic optimization, it is not possible to answer whether the global minimal was reached. Fig. 7.c depicts validation loss with the total number of neurons considered within networks, which shows a validation loss reduction for network forms of 4 to 11 neurons; however, the validation loss plateaued for networks with more than 11 neurons. As a result, increasing beyond 11 neurons does not produce a reduction in validation loss but increases computational workload. Fig. 8.a presents validation loss with activation functions. The medians are equal to 3.78 × 10 -6 for bent, 6.40 × 10 -6 for gelu, 8.26 × 10 -6 for softplus and 3.84 × 10 -6 for swish. Using bent or swish as activation functions allows a slight validation loss reduction compared to gelu or softplus. Fig. 8.b depicts the validation loss with hidden layers. The medians are equal to 4.34 × 10 -6 for one hidden layer and 6.37 × 10 -6 for two hidden layers. In this case, networks with one hidden layer allow reducing the validation loss, which is advantageous since neural networks with one hidden layer allow reducing network complexity and computational workload compared to networks with two hidden layers. Fig. 8.c depicts the validation loss and epochs and shows that increasing epoch beyond 40 is not necessary since the trained neural network median is equal to 1.84 × 10 -5 for 20 epochs and plateaued from 40 to 180 epochs between 5.00 × 10 -6 to 6.00 × 10 -6 . Fig. 9.a presents the validation loss with optimizers adam, nadam and radam. The validation loss medians are equal to 5.21 × 10 -6 for adam, 7.37 × 10 -6 for nadam and 3.50 × 10 -6 for radam. The radam reduced median validation loss compared to adam and nadam; however, because nadam and radam are derivative from the adam optimizer, the median validation loss remains similar, and no breakthroughs are noticeable. In addition, Fig. 9.b shows that the learning rate lower than 2.50 × 10 -3 allows reducing the median validation loss compared to learning rates higher than 2.50 × 10 -3 . Furthermore, Fig. 9.c and Fig. 9.d show that momentums exceeding 0.9 allow reducing the validation loss. The median is equal to 2.49 × 10 -6 for the momentum exponential decay compared to 1.53 × 10 -5 and 6.71 × 10 -6 . Additionally, median is equal to 3.01 × 10 -6 for momentum estimate compared to 1.35 × 10 -5 and 5.07 × 10 -6 .

Table 6 presents the top 10 with trained neural networks which have the lowest validation loss, which shows that both one and two hidden layers are represented; however, there is only one network with two hidden layers. The number of neurons goes from 11 to 36, where increasing the number of neurons is not necessary, as shown in Fig. 7.c. The fourth activation functions are present within the top 10; however, there is a majority of swish with 6 over 10. In addition, the epoch number from the top 10 confirms the information extracted from Fig. 8.c, showing that increasing the number of the epoch is not necessary since there are trained networks with 40 or 80 epochs. Furthermore, radam and nadam prevail as optimizers, and adam is not represented within the top 10. Regarding learning rate, the results from Fig. 9.b are confirmed since a majority of the top 10 were trained with learning rates lower than 2.50×10 -3 . Concerning momentums, the results from Fig. 9.c and Fig. 9.d are confirmed since most networks were trained with momentums higher than 0.9. Furthermore, the top 10 trained network validation loss goes from 9.46 × 10 -7 to 8.11 × 10 -7 , which is close. Finally, we selected the top network for the EMPC (highlighted in bold in Table 6), and we calculated the trained loss and test loss. The results show a training loss equal to 8.22 × 10 -7 and test loss equal to 8.38 × 10 -7 , and the network is not overfitted as test loss slightly increased compared to validation loss.

Process control

Fig. 10.a depicts the simulation results with h 1 and h 2 water level controlled by the EMPC. Water levels h 1 and h 2 are controlled in an economic manner according to the economic cost given in Eq. ( 8), without using any explicit reference trajectory. The results show that the tanks' water levels stay within the constraints. In addition, h 1 and h 2 remain close when κ i = (1; 500) and κ i = (1; 1000). Furthermore, when the cost function weights change from κ i = (1; 500) to κ i = (1; 1000) at 1,500s, the water levels h 1 and h 2 are both increased, which is the expected behavior. Moreover, h 1 and h 2 pull apart as expected when using distinct weights associated with the different levels κ i = (5; 100) or κ i = (0.5; 500). Fig. 10.b depicts the simulation results with h 3 and h 4 water level controlled by EMPC. Once again, water levels h 3 and h 4 are controlled economically within the economic cost function defined by Eq. ( 8). The latitude of permissible water levels for these reservoirs is exploited while remaining within the established limits. It can be observed that h 3 is equal to the minimum water level constraint from 3,000 to 4,500s. In addition, h 3 and h 4 remain close when κ i = (1; 500) and κ i = (1; 1000) and h 3 and h 4 pull apart when κ i = (5; 100) and κ i = (0.5; 500). Fig. 11.a depicts the q a pump flow controlled by the EMPC. During the simulation, the pump flow stays within constraints. When the cost index κ i changes from (1;500) to (1;1000), the pump flow increases. In addition, when the cost index κ i change from (1;1000) to (5;100), the pump flow decreases. However, at 3,250s, it can be seen a peak on the pump flow q a . This behavior is related to the h 3 state's constraint, during which time the pump flow increases to avoid violating the water level constraint. Finally, when the cost index κ i changes from (5;100) to (0.5;500), pump q a flow increases to maximum flow and then decreases and plateaus until the end of the simulation. Fig. 11.b depicts the q b pump flow controlled by the EMPC. During all simulations, the pump flow stays within constraints. When the cost index κ i changes from (1;500) to (1;1000), the pump flow increases. In addition, when the cost index κ i changes from (1,1000) to (5,100), the pump flow decreases to 0 m 3 /h. At 3250s however, it can be seen a peak on the pump flow q b . This behavior, identical to the previous case, is due to the state constraint h 3 , during which time the pump flow increases to avoid breaking the water level constraint. Finally, when the cost index κ i changes from (5;100) to (0.5;500), pump q b flow increases to maximum flow and then decreases and plateaus at 2.4 m 3 /h until the end of the simulation, compared to q a which plateaus at 1 m 3 /h.

Conclusion

This study aimed to propose an integrated methodology for non-linear predictive control design, synthesized from a neural network based model identified thanks to simulated data from multi-physics modelling. Results showed that the Quadruple Tank Process benchmark considered was well identified thanks to a neural network with optimized hyperparameters. Simulation using the predictive controller designed thanks to the methodology proposed depicts the expected behavior, in total consistency with the economic cost used. The presented methodology features interesting characteristics, such as being based on a fundamental movement in modeling, which consists of using bricks from existing physical models to constitute the model of the process considered. This supposes that such bricks exist, but this assumption becomes more reasonable every day, since a large community enriches the existing library.Based on this observation, the goal is to reduce the experiments providing data likely to calibrate the model, which thus performs appropriate simulations to produce the ad hoc data. On this basis, the model identifies an FNN which approaches the behavior of the physical process as closely as possible, while being less expensive to simulate. Choosing the characteristics of the target neural network (architecture and components among others) is effectively achieved by proceeding to the automated choice of the hyperparameters of the neural network by means of an adapted heuristic. Finally, the desired predictive control is obtained by choosing an economic criterion (EMPC), which is efficiently optimized thanks to choosing an ad hoc optimizer and using the FNN as an internal model. This method has been successfully implemented to control the QTP benchmark, for which different MPC solutions have already been tested in the control literature.

Among the perspectives, we are currently working to demonstrate the interest of this methodology applied to the generic problem of controlling multienergy network systems [START_REF] O'malley | Multicarrier energy systems: shaping our energy future[END_REF]. One perspective in order to avoid possible risks for the MPC control, i) to be satisfied with local minima and ii) to accelerate the resolution consists, in reformulating the optimization problem in convex and conic constraints with a disciplined convex programming [START_REF] Udell | Convex optimization in Julia[END_REF]; in this case, restrictions must be added to the neural network [START_REF] Chen | Optimal Control Via Neural Networks: A Convex Approach[END_REF].

Figure 1 :

 1 Figure 1: Schematic view of the methodology.

Figure 2 :

 2 Figure 2: System input-output representation with inputs ū(t) and outputs x(t).

Figure 4 :

 4 Figure 4: Four tank benchmark schematic view.

Figure 5 :

 5 Figure 5: Graphical representation of the four-tank benchmark within Modelica in Dymola Software.

Figure 6 :

 6 Figure 6: Controller and plant control in simulation.

Figure 7 :

 7 Figure 7: Cost function dependency on (a) the trained networks, (b) the trained networks with clusters of 10, (c) the neuron number.

Figure 8 :

 8 Figure 8: Cost function dependency on (a) the activation function considered, (b) the hidden layers number, (c) the epoch number.

Figure 9 :

 9 Figure 9: Cost function dependency on (a) the optimizer considered, (b) the learning rate, (c) the momentum exponential decay, (d) the momentum estimate.

Figure 10 :

 10 Figure 10: Simulation results. (a) Depicts tank 1 and 2 water levels over time. (b) Depicts tank 3 and 4 water levels over time.

Figure 11 :

 11 Figure 11: Simulation results. (a) Depicts pump a flow control over time. (b) Depicts pump b flow control over time.

Table 1 :

 1 QTP parameters[START_REF] Alvarado | A comparative analysis of distributed MPC techniques applied to the HD-MPC four-tank benchmark[END_REF].

	Parameters Value		Units Description
	h 1max	1.36		m	Maximum water level tank 1
	h 2max	1.36		m	Maximum water level tank 2
	h 3max	1.30		m	Maximum water level tank 3
	h 4max	1.30		m	Maximum water level tank 4
	h min	0.2		m	Minimum water level all tanks
	q amax	3.26		m 3 /h Maximum water flow pump a
	q bmax	4.00		m 3 /h Maximum water flow pump b
	q min	0		m 3 /h Minimum water flow all pumps
	S c	0.06		m	All tanks cross-section
	a 1	1.31	×	m 2	Surface of the tank 1 leakage orifice
		10 -4
	a 2	1.51	×	m 2	Surface of the tank 2 leakage orifice
		10 -4
	a 3	9.27	×	m 2	Surface of the tank 3 leakage orifice
		10 -5
	a 4	8.82	×	m 2	Surface of the tank 4 leakage orifice
		10 -5
	V a	0.3		-	Three way valve a opening
	V b	0.4		-	Three way valve b opening

Table 2 :

 2 Parameters of the PRBS signals.

	Parameters	Value	Units
	Time period	500	s
	Start time	500	s
	End time	30	days
	Amplitude

Table 3 :

 3 Parameters of constant input signals.

	Parameters	Value	Units
	Time period	1200	s
	Start time	30	days
	End time	600	days
	Amplitude

Table 4 :

 4 Hyperparameters range values.

	Symbol	Description	Range value
	σ	Activation function	σ ∈ [Bent [59], Gelu [28], Soft-
			plus [23], Swish [54]]
	δ	Neurons number	δ ∈ [4 to 19]
	θ	Hidden layer number	θ ∈ [1,2]
	µ	Epoch number	µ ∈ [20 to 180]
	ν	Batch size	ν ∈ [2048]
	o	Optimizer	o ∈ [adam [38], radam [44],
			nadam [15]]
	τ	Learning rate	τ ∈ [1 × 10 -7 to 0.01]
	υ	Momentum exponential decay υ ∈ [0.7 to 0.999]
	φ	Momentum estimate	φ ∈ [0.7 to 0.999]

Table 5 :

 5 κ and optimal steady state values.

	Time	κ i	x * 1 ; x * 2 ; x * 3 ; x * 4
	1 to 1,500	(1; 500) (0.519; 0.524; 0.514; 0.543)
	1,500 to 3,000 (1; 1000) (0.689; 0.675; 0.705; 0.668)
	3,000 to 4,500 (5; 100) (0.244; 0.291; 0.203; 0.370)
	4,500 to 6,000 (0.5; 500) (0.746; 0.512; 1.031; 0.260)

Table 6 :

 6 Hyperparameters optimization results top 10 validation loss.

	-	Hidden	Amount	Activation	Epoch	Optimizer	Learning	Momentum	Momentum	Validation
		layers	of neu-	function			rate	exponential	estimate	loss
			rons					decay
	1	1	19	bent	180	nadam	2.18 × 10 -3 9.97 × 10 -1 8.97 × 10 -1 8.11 × 10 -7
	2	1	17	swish	140	radam	2.71 ×10 -3	9.79 ×10 -1	9.22 ×10 -1	8.22 × 10 -7
	3	1	12	bent	180	nadam	1.12 ×10 -4	9.82 ×10 -1	8.52 ×10 -1	8.93 × 10 -7
	4	1	16	swish	100	radam	3.08 ×10 -3	9.79 ×10 -1	9.22 ×10 -1	9.14 × 10 -7
	5	1	11	swish	140	nadam	2.71 ×10 -3	9.79 ×10 -1	9.22 ×10 -1	9.18 × 10 -7
	6	1	11	swish	140	nadam	2.71 ×10 -3	9.79 ×10 -1	9.22 ×10 -1	9.19 × 10 -7
	7	2	36	swish	40	radam	4.55 ×10 -4	9.72 ×10 -1	9.09 ×10 -1	9.35 × 10 -7
	8	1	17	swish	80	radam	3.74 ×10 -3	9.84 ×10 -1	9.72 ×10 -1	9.35 × 10 -7
	9	1	14	softplus	160	radam	5.58 ×10 -3	9.73 ×10 -1	9.82 ×10 -1	9.43 × 10 -7
	10	1	16	gelu	180	radam	3.16 ×10 -4	9.07 ×10 -1	9.32 ×10 -1	9.46 × 10 -7

Email addresses: pierre-clement.blaud@imt-atlantique.fr (Pierre Clément Blaud), philippe.chevrel@imt-atlantique.fr (Philippe Chevrel) Preprint submitted to ElsevierDecember 12, 2021

Acknowledgment

This research was funded by Conseil Régional des Pays de la Loire, France. The authors thank the Julia community for online comments about the Julia language and packages.

Data availability statement

Research data are not shared.


Ghina Dandachi

Yassine Hadjadj-Aoul

Abdelkader Outtagarts

On Combining Reinforcement Learning and Monte Carlo for Dynamic Virtual Network Embedding

Keywords: Virtual network embedding, Deep Q-Network, Control theory

Network slicing is one of the key building blocks in the evolution towards "zero touch networks". Indeed, this will allow 5G and beyond 5G networks to deploy services dynamically, on the same substrate network, regardless of their constraints. In this demo, we introduced a platform for dynamic virtual network embedding, a problem class known to be NP-hard. The proposed solution is based on a combination of a deep reinforcement learning strategy and a Monte Carlo (MC) approach. The idea here is to learn to generate, using a Deep Q-Network (DQN), a distribution of the placement solution, on which a MC-based search technique is applied. This makes the agent's exploration of the solution space more efficient.

I. INTRODUCTION

Networks' operators (NOs) have been focusing for several years on evolving their networks to accommodate the current and future services, with reducing operational and investment costs. Among these evolutions, the virtualization of network functions (VNFs) [START_REF] Yi | A comprehensive survey of network function virtualization[END_REF], and more generally services, is the one that has brought the most in-depth changes for network management. The concept of Network Slicing (NS) has been introduced in 5G to support a wide range of services within the same substrate network.

The problem of network slicing refers to the placement of constrained services [START_REF] Rost | Network slicing to enable scalability and flexibility in 5g mobile networks[END_REF]. This research problem generalizes the bin-packing problem, which has been extensively studied in the literature and well known to be an NP-hard problem. Indeed, service placement, as envisioned by 5G, implies not only the placement of nodes, representing the components of the service, but also links representing the network constraints existing between these different components (i.e., services in the form of a graph, namely Virtual Network Function Forwarding Graph (VNF-FG)).

Several approaches have been proposed in the literature to solve the service placement problem [START_REF] Redana | 5G PPP Architecture Working Group -View on 5G Architecture[END_REF]. In [START_REF] Quang | A deep reinforcement learning approach for vnf forwarding graph embedding[END_REF], the authors proposed to combine a Deep Reinforcement Learning (DRL) strategy with a heuristic, to make the placement safer at the cost, however, of effectiveness. In [START_REF] Haeri | Virtual network embedding via monte carlo tree search[END_REF], a Monte Carlo Tree Search (MCTS) strategy is proposed. It allows to find a suboptimal solution to the placement problem, whereas the cost of a new research remains substantial since there is no learning. More recently, the authors, in [START_REF] Yan | Automatic virtual network embedding: A deep reinforcement learning approach with graph convolutional networks[END_REF], proposed to combine Graph Neural Networks with DRL to efficiently solve the service embedding problem.

Unlike the majority of existing works, we present, in this demo paper, a demonstrator of a DRL technique that considers dynamic requests of network slices with dynamic requirements and resource requests. The proposed solution leverages on DQN [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF] with Monte Carlo (MC) search to learn the best placement solutions. The combination of these two approaches increases the exploration of the solutions' space. The proposed demo shows that our proposal improves the performance of heuristics such as a First-Fit, but also to improve the approach based on a pure MC strategy.

The rest of this paper is structured as follows. In Section II, we introduce the architecture of the proposed solution. On Section III, we introduce the experimentation settings and the obtained results. Section IV concludes the paper.

II. SYSTEM ARCHITECTURE

The proposed simulator demo is implemented in Python, using the PyTorch library. Figure 2 shows the GUI of the demonstrator, the user can choose the simulation time, the strategy, and the number of MC iterations. The problem resides in the ability of the controller to learn while having dynamic arrival of slices' requests, contrarily to most existing work, considering a fixed setting for slices' requests.

A. Architecture description

The architecture of our demo is presented in Figure 1. The Simulator loads the substrate information and receives dynamic virtual network requests (VNRs). The DRL agent extracts the features from the system state defined by the substrate available resources and VNR requirements. The MC strategy allows the explorer to generate several node-by-node DQN-based distributions for the VNR placement. The explorer selects a feasible placement with the highest revenue to the cost (r2c) and checks the nodes placement feasibility as well as edges connecting the VNR nodes. Then, it updates the substrate information with the new VNR placement. The DQN agent writes the placement solution to the batch memory for future learning.

1) System State:

• Substrate state: Defined in the file Substratesim.py to create the substrate instance and update its available resources.

• VNR state: Defined in vnrsim.py to create an instance of an arriving VNR with its requirements. extracts substrate nodes selection probability distribution using DQN selects a substrate node for the current VNR node obtain a VNF-FG and calculates its corresponding r2c value keeps track of all the states visited during the selected path for eventual usage by the DQN agent if the path is chosen to be the best path • selects the VNF-FG with the highest r2c 4) Parameters: The Simulator parameters are defined in the file parameters.json and used during runtime. The main constants are:

• The simulation time, i.e, the number of episodes (SIM TIME).

• The explorer strategy: First-Fit (FF), MC, or DQN (implicitly with MC if the number of iterations N iter > 1). • The number of MC iterations.

• The mean time between arrival for two VNR requests (MTBA). • The minimum and maximum requested CPU (vcpu range). • The minimum and maximum request bandwidth (vbw range).

• The maximum number of VNFs in a VNR (vnfs range).

• The minimum and maximum VNR duration (vnr duration min, vnr duration max).

III. EXPERIMENTS

All the modules of our demonstrator have been installed on a laptop. We consider for this demo the placement of VNRs using the following strategies: the First-Fit heuristic, the MC heuristic algorithm and the proposed DQN with MC solution. The First-Fit heuristic consists in placing the virtual nodes in the first available physical node, and then placing the paths using the Dijkstra algorithm.

To assess the performance of the proposed approach, we consider a network topology with 5 nodes and 4 full-duplex links. We considered the placement of dynamic requests of VNF-FGs, comprising between 3 VNFs each. For the physical nodes, we considered CPU in the interval [20,50]. For physical links, we have considered bandwidth in the interval [20,50]. For the virtual nodes, we considered CPU resources in the interval [START_REF] Haeri | Virtual network embedding via monte carlo tree search[END_REF]10]. For the virtual links we considered bandwidth in the interval [START_REF] Haeri | Virtual network embedding via monte carlo tree search[END_REF]10].

Figure 3 shows the results obtained with different strategies: First-Fit (FF), DQN without MC (DQN), MC, and DQN with MC (DQN + MC) for different number of iterations (2 and 4 iterations). The obtained results show the revenue to the cost observed by the controller for a VNR placement defined as the amount of resources requested by a VNR divided by the actual amount of resources needed to achieve a VNR placement. We observe that DQN without MC achieves higher r2c than FF. By adding MC iterations to the path selection, the controller is able to reach higher values of r2c.

IV. CONCLUSION

Service placement is an important technology required by 5G networks, especially for Intent based orchestration that considers dynamic patterns as well as autonomous service placement and control. In this demo, we propose a dynamic slice placement tool for network slices based on deep reinforcement learning. It builds over the previously placed slicing, and offers a long term gain when compared to the baseline strategy (First-Fit). It also guarantees a minimum performance similar to that offered by First-Fit. In the future, we are planning to train our agent on a larger substrate network topology and larger VNRs. We also intend to extend our model definition and training algorithm to support the energy status of substrate nodes in order to detect nodes failure during embedding.

Fig. 1 .

 1 Fig. 1. Demo architecture

Fig. 2 Fig. 3 .

 23 Fig. 2. Demo GUI


Panagiotis Sapountzis

Serafino Teseo

Saria Otani

Frank Møller Aarestrup

Evelyne Forano

Garett Suen

George Tsiamis

Bradd Haley

Jo Ann Van Kessel

Sharon A Huws

The Fecobiome Initiative

Keywords: research initiative, gastrointestinal tract, microbiota, zoonoses

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Background

C attle have been our domesticated associates since the dawn of agricultural societies, ca 10,000 years ago. With the complementary contributions of a highly specialized gut microbiota, they efficiently convert plant material into protein mass, utilizing dietary components that cannot be digested by monogastrics. This ability has made them our valuable partners since their domestication, as they provide us with a stable supply of protein in the form of meat and dairy products.

As expected, the microbial community residing in the gastrointestinal tract (GIT) of cattle has been the focus of research efforts for more than half a century. The rumen is the primary fermentation site (Fig. 1) where a complex microbial community, composed of archaea, bacteria, protozoa, and fungi, is responsible for the degradation and fermentation of dietary components and the provision of digested nutrients to the host [START_REF] Koike | Fibrolytic rumen bacteria: Their ecology and functions[END_REF][START_REF] Dodd | Xylan degradation, a metabolic property shared by rumen and human colonic Bacteroidetes[END_REF][START_REF] Comtet-Marre | Metatranscriptomics reveals the active bacterial and eukaryotic fibrolytic communities in the rumen of dairy cow fed a mixed diet[END_REF][START_REF] Huws | Addressing global ruminant agricultural challenges through understanding the rumen microbiome: Past, present, and future[END_REF][START_REF] Seshadri | Cultivation and sequencing of rumen microbiome members from the Hungate1000 Collection[END_REF][START_REF] Stewart | Assembly of 913 microbial genomes from metagenomic sequencing of the cow rumen[END_REF][START_REF] Kameshwar | Metadata analysis approaches for understanding and improving the functional involvement of rumen microbial consortium in digestion and metabolism of plant biomass[END_REF][START_REF] Mizrahi | The rumen microbiome: Balancing food security and environmental impacts[END_REF].

When partially digested feed and microbial biomass pass through the rumen, digestion continues in the other GIT compartments, until the contents reach the large intestine (Fig. 1), the last section of the GIT. Here, water and nutrients are absorbed [START_REF] Moran | How the rumen works. Tropical Dairy Farming Feeding Management for Small Holder Dairy Farmers in the Humid Tropics[END_REF][START_REF] Bergmann | Microbial community composition along the digestive tract in forage-and grain-fed bison[END_REF], whereas a prominent microbial biomass (second only to that of the foregut) further ferments the remaining dietary components.

The large intestine is of crucial interest not only because of the metabolic processes in which its resident microbiota engage, but also because it is the primary colonization site of zoonotic and animal pathogens. Cattle, for example, are the primary reservoir of the Shiga toxin-producing Escherichia coli (STEC). Because STEC pose no threat to an adult host, they do not trigger any host defenses; they exploit the commensal microbiota to transit through the rumen and subsequently colonize the rectal mucosa of the large intestine [START_REF] Naylor | Lymphoid follicle-dense mucosa at the terminal rectum is the principal site of colonization of enterohemorrhagic Escherichia coli O157:H7 in the bovine host[END_REF][START_REF] Sapountzis | An overview of the elusive passenger in the gastrointestinal tract of cattle: The Shiga toxin-producing Escherichia coli[END_REF]. Several 16S rRNA amplicon sequencing studies have suggested that members of the cattle rectal microbiota may be able to exclude zoonotic STEC based on mutual exclusion patterns between the rectal microbial community and STEC [START_REF] Xu | Escherichia coli O157:H7 super-shedder and nonshedder feedlot steers harbour distinct fecal bacterial communities[END_REF][START_REF] Zaheer | Changes in bacterial community composition of Escherichia coli O157:H7 super-shedder cattle occur in the lower intestine[END_REF][START_REF] Stenkamp-Strahm | Associations between Escherichia coli O157 shedding and the faecal microbiota of dairy cows[END_REF][START_REF] Wang | Interactions of the hindgut mucosa-associated microbiome with its host regulate shedding of Escherichia coli O157:H7 by cattle[END_REF]. Other E. coli zoonotic pathotypes include the enterotoxigenic E. coli, which colonize the small intestine (Fig. 1) and can be detected in animal feces [START_REF] Nagy | Enterotoxigenic Escherichia coli (ETEC) in farm animals[END_REF].

Similar to Escherichia, Salmonella pathotypes infect cattle through the fecal-oral route [START_REF] Holschbach | Salmonella in dairy cattle[END_REF] and colonize the rectal mucosa [START_REF] Palmer | Comparison of rectal mucosal cultures and fecal cultures in detecting Salmonella infection in horses and cattle[END_REF]; the establishment of several pathotypes has been suggested to be prevented by a healthy cattle GIT microbiota [START_REF] Holschbach | Salmonella in dairy cattle[END_REF], although experimental evidence is lacking. Other known zoonotics of ruminants include Clostridium perfringens [START_REF] Fohler | Diversity of Clostridium perfringens toxingenotypes from dairy farms[END_REF], Listeria monocytogenes, and Campylobacter species [START_REF] Hannon | Prevalence and risk factor investigation of Campylobacter species in beef cattle feces from seven large commercial feedlots in Alberta, Canada[END_REF].

In addition to zoonotic pathogens, which are often but not always, harmless hitchhikers of the bovine host, the large intestine is also an Achilles heel for young bovines, which are susceptible to bacterial infections. Early in life, the epimural microbial communities of the intestine are critical for the development and homeostasis of the GIT immune system [START_REF] Malmuthuge | Understanding the gut microbiome of dairy calves: Opportunities to improve early-life gut health[END_REF][START_REF] Huws | Addressing global ruminant agricultural challenges through understanding the rumen microbiome: Past, present, and future[END_REF], whereas certain members, such as lactic acid bacteria, which are dominant in very young calves [START_REF] Alipour | The composition of the perinatal intestinal microbiota in cattle[END_REF], appear to protect against enteric infections and gut dysbiosis [START_REF] Malmuthuge | Understanding the gut microbiome of dairy calves: Opportunities to improve early-life gut health[END_REF][START_REF] Morgavi | Colonization and establishment of the rumen microbiota-Opportunities to influence productivity and methane emissions[END_REF].

The rapidly developing field of genomics has played an important role in deciphering functions of the symbiotic GIT microbes of animals, including humans, as well as the monitoring of antimicrobial resistance (AMR) and virulence genes. Characterization of the genetic machinery of gut microbes has FIG. 1. Overview of the FI vision and goals, the relevant research areas, and how successful implementation could improve animal health and nutrition and reduce the zoonotic and AMR exposure risk to humans. Black arrows show benefits. On top (highlighted with a light blue background), the two main microbial groups of interest (separated based on whether they are beneficial mutualists or animal/zoonotic pathogens) and their implications to potential human exposure and animal husbandry practices are presented. AMR, antimicrobial resistance; FI, Fecobiome Initiative.

improved our understanding of: (i) mutualistic (typically nutritional) services gut microbes provide to their hosts; (ii) microbial interactions; and (iii) proximate and evolutionary mechanisms shaping genomes, including factors that facilitate the development of pathogenic traits that can exploit the host and/or associated microbial communities.

Besides the nutritional aspects that may apply to almost all animals [START_REF] Bourtzis | Insect Symbiosis, 1st ed[END_REF][START_REF] Matthews | The rumen microbiome: A crucial consideration when optimising milk and meat production and nitrogen utilisation efficiency[END_REF][START_REF] Gardiner | Impact of intestinal microbiota on growth and feed efficiency in pigs: A review[END_REF][START_REF] Jahnes | Nutritional symbiosis and ecology of host-gut microbe systems in the Blattodea[END_REF], the interplay between microbiota and host has also been proven to impact all aspects of the host's health, development, hormone regulation, and behavior, as almost every physiological process is potentially affected by the microbes residing in the GIT [START_REF] Hara | The role of the intestinal microbiota in type 1 diabetes[END_REF][START_REF] Kamada | Control of pathogens and pathobionts by the gut microbiota[END_REF]Kanji et al., 2017;[START_REF] Malmuthuge | Understanding the gut microbiome of dairy calves: Opportunities to improve early-life gut health[END_REF][START_REF] Butt | Gut microbiota and energy homeostasis in fish[END_REF][START_REF] Teseo | The Scent of Symbiosis: Gut bacteria may affect social interactions in leaf-cutting ants[END_REF][START_REF] Liberti | The gut microbiota-Brain axis of insects[END_REF][START_REF] Sapountzis | An overview of the elusive passenger in the gastrointestinal tract of cattle: The Shiga toxin-producing Escherichia coli[END_REF].

For the past 10 years, the use of metagenomics and nextgeneration sequencing studies has led to important insights in all the earlier mentioned research areas, allowing researchers to: (i) take an integrated approach to the characterization of microbial gene inventories (shotgun metagenomics) and the reconstruction of metagenome-assembled genomes of uncultured symbionts, also known as ''dark matter'' [START_REF] Bernard | Microbial dark matter investigations: How microbial studies transform biological knowledge and empirically sketch a logic of scientific discovery[END_REF]; (ii) perform quantitative comparisons of the diversity and structure of microbial communities (biomarkers such as 16S, 18S, and ITS).

In cattle, genomic approaches (both culture and nonculturebased) have mostly focused on the rumen [START_REF] Pitta | Metagenomic assessment of the functional potential of the rumen microbiome in Holstein dairy cows[END_REF][START_REF] Seshadri | Cultivation and sequencing of rumen microbiome members from the Hungate1000 Collection[END_REF][START_REF] Stewart | Assembly of 913 microbial genomes from metagenomic sequencing of the cow rumen[END_REF][START_REF] Bickhart | Assignment of virus and antimicrobial resistance genes to microbial hosts in a complex microbial community by combined long-read assembly and proximity ligation[END_REF]. Despite the clear role of the large intestine in many aspects of health, nutrition, and the spread of zoonotic pathogens and AMR, the knowledge of its microbiota is mainly limited to:

(i) descriptive 16S/18S rRNA targeted sequencing (Holman and Gzyl, 2019), (ii) single bacterial species (i.e., Escherichia and Salmonella) and their antagonists that pose a threat to the cattle host or human consumers [START_REF] Wang | Use of bacteriophages to control Escherichia coli O157:H7 in domestic ruminants, meat products, and fruits and vegetables[END_REF][START_REF] Holschbach | Salmonella in dairy cattle[END_REF], and

(iii) the occurrence and distribution of AMR genes in the feces (Noelle Noyes et al., 2016a, b;[START_REF] Doster | Investigating effects of tulathromycin metaphylaxis on the fecal resistome and microbiome of commercial feedlot cattle early in the feeding period[END_REF][START_REF] Zaheer | Impact of sequencing depth on the characterization of the microbiome and resistome[END_REF][START_REF] Zaheer | Comparative diversity of microbiomes and Re-sistomes in beef feedlots, downstream environments and urban sewage influent[END_REF][START_REF] Keijser | Dose-dependent impact of oxytetracycline on the veal calf microbiome and resistome[END_REF][START_REF] Liu | The fecal resistome of dairy cattle is associated with diet during nursing[END_REF][START_REF] Rovira | Characterization of the microbial resistome in conventional and ''Raised Without Antibiotics'' beef and dairy production systems[END_REF][START_REF] Weinroth | Investigation of tylosin in feed of feedlot cattle and effects on liver abscess prevalence, and fecal and soil microbiomes and resistomes[END_REF][START_REF] Gaeta | A cross-sectional study of dairy cattle metagenomes reveals increased antimicrobial resistance in animals farmed in a heavy metal contaminated environment[END_REF][START_REF] Haley | Differences in the microbial community and resistome structures of feces from preweaned calves and lactating dairy cows in commercial dairy herds[END_REF][START_REF] Lim | Antibiotic resistomes discovered in the gut microbiomes of Korean swine and cattle[END_REF][START_REF] Salaheen | Metagenomic analysis of the microbial communities and resistomes of veal calf feces[END_REF][START_REF] Duarte | Metagenomics-based approach to source-attribution of antimicrobial resistance determinants-Identification of reservoir resistome signatures[END_REF] (see also https://www .fecobiome.com/resources/data/).

As of now, very few functional studies or shotgun metagenomic studies have taken an integrated approach to explore the genomic potential of the beneficial microbial residents, the characterization of animal pathotypes or zoonotic hitchhikers and their interactions with the host or each other, the ecological factors impacting their biology, or the potential relationship between microbiota characteristics and the spread of AMR and virulence factors.

We believe it is timely to establish an initiative that can identify common ground between these types of projects, help exploit the full potential of existing sequencing data and other resources, and provide a platform for the emergence of ideas that can generate novel functional interdisciplinary projects. Examples of such projects include the genomic and functional deciphering of the microbial di-versity of the large intestine, which can provide the cornerstone for all kinds of follow-up transcriptional, chemical, or metabolomic studies and a global effort to collect microbial isolates representing members from all genera that can be found in cattle GIT.

By helping the emergence of such projects (via experimental and intellectual collaboration), which, in turn, can lead to a better understanding of the global microbial diversity of the GIT microbiome, we are confident that the Fecobiome Initiative (FI) can ultimately limit zoonotic pathogens from entering our food chain, better control animal pathogens that hinder the animal development and productivity, reduce the emergence of AMR, optimize the use of antibiotics in animal husbandry, and potentially maximize the nutritional benefits provided by the microbial community to the cattle host.

FI Vision

We hereby propose the formation of the FI (similar to the Greek letter F; also known as the golden ratio in mathematics), an international effort that will promote scientific cooperation, dissemination of research findings, and data transparency to facilitate the advancement of research related to microbial communities primarily of the large intestine of cattle, but also related to the rest of the GIT. The vision of the FI is to offer viable solutions in cattle husbandry that will help reduce the zoonotic and AMR risk, animal pathogens, and reduce or optimize the use of antibiotics, with the goal of improving animal welfare while enhancing meat and dairy production.

FI Mission

Through a dedicated website (www.fecobiome.com) that will serve as a coordination and dissemination agent, the FI will:

(i) bring together and coordinate collaboration among members of the scientific community (e.g., by organizing workshops to align strategies for funding applications and joint proposals) who are interested in the study of microbes residing in the GIT (primarily the large intestine) of cattle, (ii) disseminate knowledge by using language and material targeting not only scientists (working on cattle microbiota-related research) but also nonspecialists (e.g., again by organizing events like webinars or workshops for training purposes) and stakeholders involved in animal husbandry, and thus help bridge the gap between the scientific community and nonsubject matter specialists, (iii) increase public availability of protocols, relevant material (e.g., primers, probes, antibodies), and isolated microbial strains; and by making data (genomic, metagenomic or other) more accessible to the public (i.e., by contributing whenever possible with supplementary metadata and a built-in interface that will allow easier search and recovery of genomic data).

Why Feces?

The organizers of the FI are interested in the study of gut microbes residing primarily in the last compartment of the cattle GIT, the large intestine. To achieve this, the FI will focus on samples originating from fresh feces as: (i) they present a reliable representation of the microbial community at the distal end of the GIT, and (ii) they do not need special permits for handling of animals in commercial or other farms. In addition, feces can also contain microbial isolates from representatives of other parts of the GIT such as the small intestine and the rumen, which will be valuable for research related to these parts of the GIT (see also the RMG network; https://rmgnetwork.org/).

Fresh feces is an excellent sample for: (i) the isolation of cultivable microbes; (ii) isolation and characterization of nucleotide or protein material for -omics approaches, diagnostic polymerase chain reaction, or other molecular and biochemical methods; (iii) chemical analyses for determining nondigested feed components; and (iv) in vivo and in vitro experimentation. We would note that, despite our interest in the microbial communities of fecal samples, FI does not argue against the use of material from the GIT of animals that has been collected via direct rectal sampling (e.g., rectal swabs) or necropsies (i.e., GIT parts from slaughtered calves) as this material will remain invaluable for other research questions.

FI Goals

The FI has the potential to benefit research in the areas of ecology, evolution, epidemiology, microbiology, human health, and medicine and veterinary medicine, and it may lead to solutions for the improvement of animal husbandry practices. Some examples include: Characterization of the distal end microbiome will allow us to formulate strong hypotheses about the nature of the relationship between members of the microbial community and the cattle host, the intertwined dynamics among the community members (cooperative or competitive interactions), and the functional complementarity or redundancy across compartments (e.g., when examined in comparison with similar available datasets from the rumen; Hungate dataset). Insights gleaned from such studies can significantly improve our understanding of the homeostasis and functioning of the GIT microbial community. Characterizing the microbial variability on a global scale will allow us to link the variability to host and environmental factors (through correlations of big data sets) such as climate, geographic location, diets, and farm management practices that may affect microbial establishment and resilience. This will allow us to formulate hypotheses about the key factors that influence the GIT microbial community. The Whole Genome Sequencing approaches may reveal patterns of genomic evolution, since comparisons of the isolated bacterial strains with their close relatives can reveal traits that have evolved as a result of the unique ecosystem in the large intestine.

The FI can directly contribute to applied solutions in animal husbandry, including the discovery of novel direct-fed microbials to potentially protect against pathogens during early development, to prevent zoonotic pathogens from entering and establishing themselves in the GIT, to limit the spread of AMR bacteria, and potentially to improve animal nutrition and productivity.

In addition, the FI could lead to improvements in the assessment of zoonotic risks, since currently zoonotic risk and the presence of pathogens have been evaluated by targeting only a handful of known pathogens (mainly Escherichia and Salmonella genera) and with culture-based approaches. By promoting culture-independent approaches that examine the entire fecal community in parallel with culture-based approaches, we will enhance our ability to identify emerging pathogens or other threats.

Short-Term and Long-Term Objectives and Benefits to the Community

The FI will identify knowledge gaps and focus on research that will complement existing knowledge regarding the microbial inhabitants of the large intestine.

To achieve this, the FI will prioritize research areas such as:

(i) Characterization of the microbial community of cattle feces by using high-quality/in-depth shotgun metagenomic sequencing to make detailed gene catalogs of the microbial genomic inventory present in the large intestine and identify potential pathogenic and AMRassociated genes. (ii) Isolation of microbial strains from fecal samples using both taxa-specific and nonspecific approaches and create a representative library of native aerobic and anaerobic isolates. (iii) Using complementary -omics approaches to assemble evidence-based data as these approaches can link, among other factors, DNA sequences, proteins, and metabolites. (iv) Large-scale correlations between (meta) genomic data to identify potential factors responsible for the variations between microbial communities (e.g., ecological and geographic datasets, diets, husbandry practices, breed, age, sex, etc.). These correlations require data and isolates from cattle of every breed and age around the globe. We recognize this is an ambitious goal, but we are hopeful that, once the FI is launched, the FI can facilitate data collection, leading to novel results that fuel future projects that allow us to achieve our goals.

We note that the dissemination of data and knowledge is one of the core missions of the FI and toward this end we plan to create a wiki portal (as part of the fecobiome.com website) that will present relevant scientific articles in a simplified and concise view to allow nonsubject matter experts to better connect with the scientific community. Bridging the gap between animal husbandry stakeholders and researchers has the potential to significantly increase the availability of fecal samples from farms, as currently the communication gap between these two parties is often an obstacle.

As (meta)genomic data accumulate, the map of fecal microbial diversity will become more complete, and accomplishing the FI short-and long-term objectives will aid all projects related to the GIT distal end microbiome: Microbial isolations and genomic sequencing will provide reference genomes that will serve as templates for meta-genomic, transcriptomic, and proteomic studies and the identification of novel key genes and taxa (e.g., to be used as markers). This includes pathogenic taxa and virulence genes that can be used as biomarkers, as global comparisons will allow for the evaluation of the most relevant threats. Complementary in vitro and in vivo methods applied to microbial isolates can evaluate their primary and secondary metabolic potential and growth dynamics; insight that, when coupled with metagenomic data, can provide a reconstruction of the community metabolic network and potential ecological interactions. Characterizing the large intestine microbes will foster synergies with existing networks related to the cattle GIT microbiome (i.e., RMG network) and ongoing projects, as the combined results can provide an overview of the community structure and mechanisms along the GIT and allow for an integrated approach to the study of the GIT.

Launching the FI

Considering the growing number of shotgun metagenomic fecal data in cattle during the past 3 years (Noelle Noyes et al., 2016a, b;[START_REF] Doster | Investigating effects of tulathromycin metaphylaxis on the fecal resistome and microbiome of commercial feedlot cattle early in the feeding period[END_REF][START_REF] Zaheer | Impact of sequencing depth on the characterization of the microbiome and resistome[END_REF][START_REF] Zaheer | Comparative diversity of microbiomes and Re-sistomes in beef feedlots, downstream environments and urban sewage influent[END_REF][START_REF] Liu | The fecal resistome of dairy cattle is associated with diet during nursing[END_REF][START_REF] Rovira | Characterization of the microbial resistome in conventional and ''Raised Without Antibiotics'' beef and dairy production systems[END_REF][START_REF] Salaheen | Metagenomic analysis of the fecal microbiomes from Escherichia coli O157:H7-shedding and nonshedding cows on a single dairy farm[END_REF][START_REF] Salaheen | Metagenomic analysis of the microbial communities and resistomes of veal calf feces[END_REF][START_REF] Weinroth | Investigation of tylosin in feed of feedlot cattle and effects on liver abscess prevalence, and fecal and soil microbiomes and resistomes[END_REF][START_REF] Gaeta | A cross-sectional study of dairy cattle metagenomes reveals increased antimicrobial resistance in animals farmed in a heavy metal contaminated environment[END_REF][START_REF] Haley | Differences in the microbial community and resistome structures of feces from preweaned calves and lactating dairy cows in commercial dairy herds[END_REF][START_REF] Lim | Antibiotic resistomes discovered in the gut microbiomes of Korean swine and cattle[END_REF][START_REF] Duarte | Metagenomics-based approach to source-attribution of antimicrobial resistance determinants-Identification of reservoir resistome signatures[END_REF], launching the FI is timely, as its website (fecobiome.com) can serve as a central repository to organize metadata and facilitate meta-analyses to fully exploit the potential of such data.

For this purpose, we have developed a search tool at fecobiome.com (a beta version is already available online at: www.fecobiome.com/resources/) that will allow users to easily locate genomic and other data of interest. Researchers can submit the metadata of their studies through the website (https://www.fecobiome.com/contact/submit-data/). Fecal microbial isolates are already available at the National Research Institute for Agriculture, Food and the Environment (France), the Institute for Global Food Security at Queen's University Belfast (Ireland), and the Danish Technological Institute (Denmark) and the metainformation of these isolates will be shortly available at our website. Similarly, a list of protocols and reagents along with research labs that can provide them will also be made publicly available through www.fecobiome.com.

The FI website will serve as a platform for connection and communication among members (an internet forum is an upcoming feature), where members can exchange ideas related to our initiative. The website administrators will be responsible for regulating the threads and for monthly news updates. An Internet forum can also facilitate communication between members, which will initially take place through the website. However, in the coming years, the FI will actively pursue funding for organizing meetings and workshops where members can interact in person.

We invite all interested parties to join the FI. It will be a dynamic site and, through the continuous feedback of members, the FI will adapt to the needs of fecobiome related research in the coming years.

Disclosure Statement

No competing financial interests exist.

SAPOUNTZIS ET AL.

SAPOUNTZIS ET AL.

SAPOUNTZIS ET AL.

Funding Information

This publication was supported in part by a start package grant from the National Research Institute for Agriculture, Food and the Environment.