Charles-Olivier Artizzu
email: artizzu@i3s.unice.fr

Guillaume Allibert
email: allibert@i3s.unice.fr

Cédric Demonceaux
email: cedric.demonceaux@u-bourgogne.fr

Deep Reinforcement Learning with Omnidirectional Images: application to UAV Navigation in Forests

Keywords: Vision for robots, Mobile robotics, Perception systems

Deep Reinforcement Learning (DRL) is highly efficient for solving complex tasks such as drone obstacle avoidance using cameras. However, these methods are often limited by the camera perception capabilities. In this paper, we demonstrate that point-goal navigation performances can be improved by using cameras with a wider Field-Of-View (FOV). To this end, we present a DRL solution based on equirectangular images and demonstrates its relevance, especially compared to its perspective version. Several visual modalities are compared: ground truth depth, RGB, and depth directly estimated from these 360°R GB images using Deep Learning methods. Next, we propose a spherical adaptation to take into account the spherical distortions of omnidirectional images in the convolutional neural networks (CNNs) used in the actor-critic network and show a significant improvement in navigation performance. Finally, we modify the perspective depth estimation network using this spherical adaptation and demonstrate a further performance improvement.

I. INTRODUCTION

Depth is the most commonly used input for drone navigation in complex, unstructured environments. This visual modality enables reliable obstacle detection and thus safe trajectories for Unmanned Aerial Vehicles (UAVs). Most of today's commercial drones use LIDARs to obtain an accurate depth estimation. However, these sensors are heavy and power-consuming. With recent accuracy improvements in monocular perspective depth estimation, RGB cameras are now being considered to replace LIDARs.

In parallel, more affordable and accurate omnidirectional cameras have recently been commercialized. The latest generation of cameras provides spherical images with up to 4K resolution and high-frequency capture. In addition, with their 360°Field-Of-View (FOV), they can capture the entire environment in a single image. Therefore, this sensor offers great potential in many robotic tasks, including drone navigation.

To solve the point-goal navigation, most state-of-the-art algorithms use a local planner based on deep reinforcement learning (DRL). This method has been proven to be effective in driving short-term navigation, and collision avoidance in many papers [START_REF] Guastella | Learning-based methods of perception and navigation for ground vehicles in unstructured environments: A review[END_REF], [START_REF] Loquercio | Learning high-speed flight in the wild[END_REF], [START_REF] Kastner | Connecting deep-reinforcement-learning-based obstacle avoidance with conventional global planners using waypoint generators[END_REF]. The DRL approach proposes to learn a navigation policy through trial-and-error experiments where the agent interacts with the environment based on its perception, and state [START_REF] He | Explainable deep reinforcement learning for uav autonomous path planning[END_REF], [START_REF] Kahn | Badgr: An autonomous selfsupervised learning-based navigation system[END_REF]. A reward is designed to promote or prevent specific behaviors and thus influence the agent's policy. This method offers excellent generalization capabilities but requires a very long learning process, usually performed in virtual environments to allow thousands of trials and exploration of failure cases.

Perception is crucial for these DRL algorithms because the agent selects its next action based only on what it observes. Therefore, this information must be as accurate as possible. In recent years, monocular depth estimation in perspective images has seen significant progress, reaching a very high accuracy [START_REF] Wofk | Fastdepth: Fast monocular depth estimation on embedded systems[END_REF], [START_REF] Ranftl | Vision transformers for dense prediction[END_REF], [START_REF] Ranftl | Towards robust monocular depth estimation: Mixing datasets for zero-shot crossdataset transfer[END_REF].

However, when applied directly to omnidirectional images, these perspective methods suffer from a domain shift. Indeed, these images present large distortions near the polar regions, which significantly impact the accuracy of networks designed for perspective images. Some specialized networks have been proposed to estimate depth in spherical images [START_REF] Zioulis | Omnidepth: Dense depth estimation for indoors spherical panoramas[END_REF], [START_REF] Feng | 360 depth estimation in the wild -the depth360 dataset and the segfuse network[END_REF]. In this paper, we offer an alternative method based on convolutions taking into account distortions [START_REF] Artizzu | Omniflownet: a perspective neural network adaptation for optical flow estimation in omnidirectional images[END_REF].

The contributions of this paper are fourthfold:

• First, we propose a point-goal navigation DRL solution based on 360°FOV images as input.

• Second, we compare the performance of this solution with its baselines using perspective images. To our knowledge, this is the first comparison of omnidirectional and conventional images for DRL. The obtained results confirm the advantage of omnidirectional systems.

• Third, the ground truth depth, RGB, and estimated depth modalities are compared in a performance benchmark. This comparison reveals that using the MIDAS network [START_REF] Ranftl | Towards robust monocular depth estimation: Mixing datasets for zero-shot crossdataset transfer[END_REF] is more reliable than methods based on RGB alone.

• Fourth, we propose an adaptation to take into account the spherical distortions of equirectangular images. We implement this modification in two different locations. One in the Actor-Critic network proves that spherical networks' performance is improved when trained with a specific spherical adaptation. The other directly in the MIDAS network demonstrates that a pre-trained network can be adapted to equirectangular images without additional training.

The structure of this paper is the following. Section II describes the goal-driven navigation solution adopted here. Next, Section III presents the FOV and modality performance benchmark. Finally, Section IV presents the proposed adaptation for spherical images and the associated results.

II. DRL-BASED NAVIGATION SOLUTION

Point-goal navigation can be modeled by a Markov decision process (MDP). In this formulation, the agent interacts with the environment by performing actions following a specific policy in a given state. The agent receives a positive or negative reward, to promote or prevent certain behaviors. In our specific case, the drone is rewarded positively when it reaches its goal and negatively when it collides with obstacles, gets stuck in a loop, or moves too far away from its objective.

A. State

In DRL, the agent uses its current state to determine its following action. Therefore, the observation of his state must be accurate and complete enough to provide sufficient information to make an appropriate decision. But in return, a too exhaustive state will overload the agent with redundant parameters. Thus, in this study, we propose to use a state containing bare minimum information to achieve the two main objectives: point-goal navigation and obstacle avoidance.

For the navigation task, only the relative distance and direction of the current goal are provided. At each time-step t k = k∆t, where ∆t is the control sampling time, we define the distance d k and the angle θ k to the goal:

d k = ∥P k -P ⋆ ∥ 2 , θ k = arctan2 (y ⋆ -y k , x ⋆ -x k ) -ψ k , (1)
with a drone at position P k = (x k , y k , z k ) heading towards a fixed goal at position P ⋆ = (x ⋆ , y ⋆ , z ⋆ ) with a yaw angle ψ k . At the same time, the drone captures its surroundings with its perception sensor and transforms it into an image noted I k . Finally the state of the drone is defined at each time-step t k by:

S k = [d k , θ k , I k ].

B. Visual modalities compared

We propose to compare three different visual modalities for the I k image: Ground Truth (GT) depth, RGB, and Estimated Depth (ED) using Deep Learning. In parallel, two different capture FOVs are compared for all these modalities: a limited one (90°) and an omnidirectional one (360°). To compare these features fairly, we keep the same network architecture between the different study cases. Thus, the image resolution I k is independent of the modality or FOV and fixed at 100 × 100 pixels. Fig. 1 shows the different inputs used in this paper.

Equirectangular Perspective

C. Action

In forest point-goal navigation, the drone must navigate between trees without collision. Due to the mostly vertical structure of the trees, we propose to realize obstacle avoidance in an iso-altitude plane. Therefore, the drone controller keeps the drone's altitude constant during its flight while the pilot focuses on rotational movements. The agent moves by selecting an action a k ∈ [-π, π] that corresponds to a desired rotation angle among N b directions, which is later sent to the low-level controller [START_REF] Shah | Airsim: High-fidelity visual and physical simulation for autonomous vehicles[END_REF]:

a k = -π + 2i N b -1 π i ∈ {0, .., N b -1}. (2)

D. Network and Reward

In this paper, we focus on the perception of the DRL algorithm using different FOVs and modalities as input. Therefore, we have based ourselves on the rest of the DRL solution from already proven contributions, particularly for the actor-critic network and reward function.

The PPO algorithm [START_REF] Schulman | Proximal policy optimization algorithms[END_REF] processes the agent state using an actor-critic network to predict the next best action. The architecture proposed in this paper is based on contributions that have already proven effective for UAV navigation [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF], [START_REF] Loquercio | Learning high-speed flight in the wild[END_REF], [START_REF] He | Explainable deep reinforcement learning for uav autonomous path planning[END_REF]. First, the visual part I of the state S is preprocessed using a succession of convolutions (CONV) and fully connected networks (FC). The resulting 32-dimensional vector is then combined with the goal information (d k , θ k ) to determine the next action a k using another fully connected network (RFC). The global network is shown in Fig. 2 and requires 60k parameters (more details in Appendix A). The reward function is crucial in DRL. It describes how the agent performed during the episode. In [START_REF] Zhang | Autonomous navigation of uav in multi-obstacle environments based on a deep reinforcement learning approach[END_REF], the authors propose a solution that shows promising results in goal-driven navigation and obstacle avoidance. Their reward function is ideally suited for our problem. However, they use depth sensors in their model to maintain a safe distance to the obstacles during the flight. In order to be more independent to the modalities where absolute depth is not necessarily known, we have removed this safe distance in the reward function.

At each time-step t k , the reward is computed using the relative goal state (d k , θ k ), a penalization term (-0.02), and an additional termination reward if the agent reached an terminal state R end . The resulting function is :

R k = -0.1d k -0.05∥θ k ∥ -0.02 + R end . (3)
Three different terminal states are possible for each episode:

• the drone collides with an obstacle, so it receives a large negative reward of R end = -5;

• the drone reaches its goal (d k < 1m), so it receives a great positive reward of R end = 5;

• the drone is stuck in a loop (more than 200 time-steps) or going too far away from its goal (d k > 100m), so it receives a small negative reward of R end = -2.

III. FIRST RESULTS

We train and compare the performances of our proposed DRL solution in a simplified forest environment RDMAP. We choose Unreal Engine [START_REF] Games | Unreal engine[END_REF] as rendering software to build this environment. Connected to AirSim [START_REF] Shah | Airsim: High-fidelity visual and physical simulation for autonomous vehicles[END_REF], an open-source robotics plugin, the simulator can provide high-fidelity image capture and a low-level controller to stabilize a drone.

The training is performed on Nvidia Tesla-V100 graphic cards and takes about 8 to 10 hours. Inference on 300 drone trajectories takes between 120 and 180 minutes.

A. Visual modalities compared

As presented in Section II-B, three different modalities are compared in this article: ground truth depth, RGB and estimated depth with Deep Learning. In parallel, for each modality, two FOV are studied: one perspective and one omnidirectional.

The ground truth depth and RGB images are directly provided by the AirSim simulator. The GT depth value is cropped to the [0..5] meters range to remain representative of the capabilities of a small UAV-mounted LIDAR.

To estimate depth from equirectangular RGB images, we can either use specialized spherical networks or standard perspective solutions with spherical adaptations. Several omnidirectional networks have been published, such as [START_REF] Zioulis | Omnidepth: Dense depth estimation for indoors spherical panoramas[END_REF], [START_REF] Feng | 360 depth estimation in the wild -the depth360 dataset and the segfuse network[END_REF]. However, they require significant computational power, and their specific training on omnidirectional datasets makes them less oriented towards drone navigation. Therefore, we choose a lightweight perspective network MIDAS [START_REF] Ranftl | Towards robust monocular depth estimation: Mixing datasets for zero-shot crossdataset transfer[END_REF]. This Convolutional Neural Network (CNN) is one of the lightest and most accurate published perspective depth estimation networks, with already proven performance in several mobile and UAV applications [START_REF] Aleotti | Real-time single image depth perception in the wild with handheld devices[END_REF], [START_REF] Polosky | Machine learning subsystem for autonomous collision avoidance on a small uas with embedded gpu[END_REF]. We directly use its pre-trained version midas v21 small to be representative of an embedded drone solution. Solutions to adapt to equirectangular images are proposed in Section IV.

B. Simplified training and testing environment RDMAP

We build a simplified forest environment named RDMAP. This 200 × 200 meters terrain is made of many vertical cylinders placed randomly and schematizing a dense forest of tree trunks. Fig. 3 shows an overview of this simplified environment.

C. Training

The proposed DRL solution is trained in RDMAP using the different visual modalities and FOVs. Each training takes 100k time-steps and uses the same schedule and parameters for a fair comparison. Appendix A shows the hyperparameters used to tune the PPO algorithm. During training, the distance between the initial position of the drone and its target is always 20 meters.

D. Metrics

Two standard metrics are used to evaluate each solution navigation results:

• the Success Rate SR (the percentage of successful runs divided by the total number of runs) to directly assess the drone's abilities to reach its goal (d k < 1m);

• the Success weighted by Path Length, as defined in [START_REF] Anderson | On evaluation of embodied navigation agents[END_REF]:

SP L = 1 N N i=1 S i ℓ i p i ,
where S i is 0 or 1 depending on the success of the episode, p i is the length of the drone's trajectory and ℓ i is the shortest distance between the initial and goal points. Thus, the closer the drone trajectory is to the shortest path, the closer the SPL is to 100%. Besides, failed tests are strongly penalized by the boolean value S i .

E. Inference in RDMAP

We test the performance of the trained DRL solution directly in the training environment RDMAP. The testing process is identical for each modality or FOV: point-goal navigation is proven on the same 600 randomly selected drone-goal pairs on the map. However, unlike training, the distances between the targets are not only 20m: we also used distances of 40m and 60m to test our solution and prove its robustness thoroughly.

As presented before in Section II-B, three visual modalities (GT depth, RGB, and Estimated Depth (ED) with MIDAS) and two FOVs (perspective and omnidirectional) are compared. The results are presented below in Table I.

Solutions using omnidirectional images consistently outperform those using perspective ones. The larger FOV improves navigation and obstacle detection tasks. Obstacles invisible in a narrow FOV are now detected in the omnidirectional images reducing the number of collisions (higher SR). Moreover, the drone has a better understanding of its environment resulting in better trajectory optimization (higher SPL). These results demonstrate the great potential of large FOV sensors for pointgoal navigation, regardless of the type of visual information captured.

As expected, ground truth depth is the best performing visual modality. In second place, estimated depth presents a slight advantage over RGB images. As a result, we conclude that MIDAS prediction accuracy, despite is lightweight architecture, is sufficient to build a navigation solution more reliable than one based on RGB alone.

IV. DRL SPHERICAL ADAPTATION

The proposed DRL-based navigation solution shows promising results using 360°FOV observations as inputs. However, all spherical projections introduce some distortions. In particular, equirectangular images show significant distortions near the polar regions.

Therefore, we present a solution to consider these distortions by maintaining the local coherence of pixels during convolution operations. This approach is applied in two different locations in our proposed solution. First, the convolution layers of the actor-critic network used in the DRL algorithm are modified. Second, the MIDAS network is adapted to improve depth prediction in equirectangular images without additional training.

A. Distortion-aware convolution

Several methods deal with spherical distortions in Convolutional Neural Networks (CNNs). A first approach consists in modifying the entire feature maps using Fourier transforms [START_REF] Cohen | Spherical cnns[END_REF] or spherical polyhedra [START_REF] Lee | Spherephd: Applying cnns on a spherical polyhedron representation of 360°images[END_REF]. A second technique directly changes the shape of the convolutional kernels. The new kernels shapes can either be learned [START_REF] Su | Kernel transformer networks for compact spherical convolution[END_REF] or specifically adapted to the equirectangular projection [START_REF] Fernandez-Labrador | Corners for layout: End-to-end layout recovery from 360 images[END_REF], [START_REF] Artizzu | Omniflownet: a perspective neural network adaptation for optical flow estimation in omnidirectional images[END_REF]. This latter method does not require new variables or additional learning, nor a complete modification of the network architecture. For these reasons, we choose to implement this distortion-aware convolution using a pre-computed offsets table.

During training and inference of the DRL solution, these offsets are directly applied to change the point position in every kernel to project the spherical image locally onto its perspective equivalent. For example, Fig. 4 shows different distortion-aware kernel shapes in function of their position in the equirectangular image.

B. Actor-Critic Network Adaptation

The four convolutional layers of the proposed actor-critic network architecture are modified. Prior to training, the offsets tables are computed based on the resolution of the observation used and the different parameters of each adapted convolutional layer. Fig. 5 shows how this additional plugin is applied on the network architecture. The proposed navigation solution is trained in the RDMAP environment in a process similar to that presented in Section III-C.

Results in RDMAP:

The equirectangular adapted DRL solution is tested in the RDMAP environment and compared to its baseline from Section III-E. Table II shows the performances of 360°FOV navigation based on ground truth or estimated depth. The baseline results are directly reused from Table I For each modality, the distortion-aware (adapted) solutions show highly better performance. Maintaining local pixel coherence during convolutions helps the Actor-Critic network to better detect obstacles in the I k input image. As a result, the drone trajectories are better optimized, with fewer collisions and faster paths. In particular, it allows the spherically adapted ED solution to outperform the non-adapted GT depth based solution, especially in trajectory optimization.

C. MIDAS network Adaptation

To further improve the results of the previous section, we propose a fast and elegant method to increase depth estimation accuracy in equirectangular images without additional training.

MIDAS is a lightweight depth estimation network trained on perspective datasets. Therefore, its accuracy suffers from spherical distortion when applied directly to equirectangular images. The previous test section showed that the DRL solution based on 360°estimated depth performs well but is still lower than ground truth depth performances.

Therefore, we propose to improve depth prediction accuracy by applying a spherical adaptation to the MIDAS network. The decoder convolutions are modified with additional precomputed offsets as explained in Section IV-A. No additional training is performed on an equirectangular dataset as we keep the training weights of the MIDAS authors on perspective datasets.

The actor-critic network spherical adaptation remains active as it has shown excellent results in the previous section. Fig. 6 shows the new global pipeline.

Results in RDMAP :

The DRL-solution using MIDAS network adapted for equirectangular images is trained and tested in the RDMAP environment. The performance of our solution with the spherical DRL and MIDAS adaptation is better than that of the DRL-only adapted solution. As a result, the MIDAS adaptation has better depth estimation accuracy than its baseline on equirectangular images.

The depth is improved locally with better obstacle detection and thus a higher success rate. In parallel, the estimation is also globally enhanced, allowing a better understanding of the complete environment of the UAV and faster trajectories (higher SPL).

Therefore, using a larger FOV and spherical adaptations allows the final version of our solution based on estimated depth to go from an initial success rate of 70% to a final version close to 90%. Moreover, the drone trajectories also gain noticeably in speed.

V. GENERALIZATION TO A PHOTOREALISTIC FOREST

ENVIRONMENT

To reduce the gap with reality, we build the RWFOREST forest environment. Using the best rendering capabilities of Unreal Engine and forest textures from its marketplace, we create a photorealistic forest with complex lighting. Unlike the simplified trunks in RDMAP, the trees here have different sizes, branches, and dense foliage. As a result, RWFOREST is a challenging environment where the captured images I k are much more complex to analyze and translate into actions. Fig. 7 shows some observations from this environment. We directly test the previously trained solutions in this new environment to challenge the robustness of our proposed solution. No additional training was performed in RWFOREST. Instead, the pre-trained models of the section III-C were directly tested on 600 objectives of 20m, 40m, and 60m. The Table IV Although a challenging test using much more complex images than those used for training, our proposed solution still performs very well. The DRL & MIDAS adapted 360°FOV estimated depth reaches almost 85% success. This proves the robustness of our proposed solution to domain change despite the increase in observational complexity.

The MIDAS network also proves its stability in predicting consistent depth estimation in various environments without additional training. Despite its small size, it is reliable enough to build a DRL navigation solution with performance close to ground truth modalities.

VI. CONCLUSION

This paper presented a point-goal navigation solution for drones in forests using monocular omnidirectional images. We have shown that this equirectangular solution outperforms its perspective reference. Indeed, thanks to a larger FOV, the drone has a better understanding of its environment, resulting in fewer collisions and better trajectory optimization.

In parallel, we have studied a DRL solution using monocular estimated depth as visual input and compared it to ground truth depth and RGB. Despite significant spherical distortions, the MIDAS network trained on perspective images provides reliable and accurate depth prediction, sufficient to outperform a DRL solution using only RGB.

Furthermore, we have presented a fast and elegant spherical adaptation that aims to maintain local pixel coherence during convolution operations without additional training. Again, the point-goal navigation task has been further improved by implementing this adaptation in the Actor-Critic network and MIDAS. Finally, even when tested in a new and more complex environment, our solutions showed high robustness to domain shift.

Our code implementation will be open-source and available on GitHub at https://github.com/COATZ/OMNI-DRL.

Fig. 1 .

 1 Fig. 1. RDMAP: (Left: RGB, Middle: depth, Right: MIDAS depth [8]).

Fig. 2 .

 2 Fig.2. Left: At each time-step t k the agent chooses an action a k based on its state S k and its policy. This interaction with the environment results in a new state S k+1 and a reward R k+1 to evaluate the previously followed policy. Right: The visual observation I k is encoded into a 32-dimensional vector using convolutional operations (CNN) and a fully connected network (FC). Then, combining this output vector and the goal information (d k , θ k ), another fully connected network (RFC) predicts the agent's next action a k .

Fig. 3 .

 3 Fig. 3. Overview of the RDMAP environment.

Fig. 4 .

 4 Fig. 4. Example of kernels with different latitude and longitude. In blue is the center of the kernel, in green the perspective kernel and in red the adapted equirectangular one. The wider distortions are near the poles.

Fig. 5 .

 5 Fig. 5. Pre-computed spherical adapted offsets are added to the four convolution layers of the Actor-Critic network during training and testing of the DRL-solution.

Fig. 6 .

 6 Fig. 6. The convolutions of the MIDAS decoder are modified without additional network training to take into account spherical distortions for depth prediction. Therefore, the final proposed DRL pipeline has two spherical adaptations: one for the MIDAS network and one for the Actor-Critic network.

Fig. 7 .

 7 Fig. 7. RWFOREST : (Left: RGB, Middle: depth, Right: MIDAS depth [8]).

  .

	RUN	SR (%)	SP L (%)
	360°FOV GT depth (baseline)	88.8	68.6
	360°FOV GT depth (DRL adapted)	94.8	78.3
	360°FOV ED (baseline)	86.0	69.0
	360°FOV ED (DRL adapted)	89.5	73.3
	TABLE II
	PERFORMANCES IN RDMAP.

  Table III shows the performances.

	RUN	SR (%)	SP L (%)
	360°FOV ED (DRL & MIDAS adapted)	89.8	74.5
	TABLE III
	PERFORMANCES IN RDMAP.

  presents the result of the best performing solutions from the previous sections: the spherical adapted ground truth and the DRL & MIDAS adapted estimated depth.

	90°FOV GT depth	81.0	69.1
	360°FOV GT depth (DRL adapted)	89.7	76.9
	360°FOV ED (DRL & MIDAS adapted)	84.7	62.1
	TABLE IV
	PERFORMANCES IN RWFOREST.

RUN

SR (%) SP L (%)

ACKNOWLEDGMENT

This work was supported by the ANR CLARA project, grant ANR-18-CE33-0004 of the French Agence Nationale de la Recherche. This work was granted access to the HPC resources of IDRIS under the allocation AD011013128 made by GENCI.


Franco Fusco

Guillaume Allibert

Olivier Kermorgant
email: olivier.kermorgant@ls2n.fr

Philippe Martinet
email: philippe.martinet@inria.fr

Investigating the Performances of Control Parameterizations for Nonlinear Model Predictive Control

published or not. The documents may come

L'archive ouverte

I. INTRODUCTION

Model Predictive Control (MPC) is probably one of the most attractive control design methodologies, as illustrated in the state of the art made by Mayne [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. The ability of MPC to handle nonlinearities, enforce constraints and provide performance/cost trade-offs are just some of the reasons for its popularity. The MPC strategy consists in computing, at each control iteration, an optimal sequence of future commands that minimizes a given cost function. Feedback is ensured since the cost explicitly depends not only on the control sequence, but also on the current measured (or estimated) states. Only the first control is applied to the system, and the whole process is repeated at the next iteration, indefinitely. Since the control problem is written in the form of a cost minimization, constraints on both the control and the state may explicitly be considered.

Historically, MPC was developed for systems with slow dynamics and required long computation times. Nonetheless, in the last two decades a lot of work has been done to deploy predictive control also on fast systems. Some of these pose additional restrictions and challenges, as the optimization problem must be solved on embedded hardware with limited computing power. One way to improve the efficiency of MPC schemes is to reduce the complexity of the optimization problem, and several methods have been proposed in order to speed up the solving process. A first one consists in approximating the original nonlinear constrained problem to find a solution as quickly as possible, even if it might be a suboptimal one. Indeed, even in case of a coarse approximation of the optimization, global performances of the closed loop controller can be satisfactory, with computation times that are sufficiently small to allow real time implementations [START_REF] Kouvaritakis | Who needs QP for linear mpc anyway?[END_REF], [START_REF]Efficient robust predictive control[END_REF]. An alternative approach proposed in the literature is to carefully exploit the particular structure of a problem to improve the efficiency of a solver [START_REF] Ferreau | An online active set strategy to overcome the limitations of explicit MPC[END_REF]. Solutions based on offline computations have also been proposed. They explore the state space offline and find an explicit solution to the problem. On-line computations are then reduced to the evaluation of a piece-wise linear function [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF].

Control parameterizations can also be an effective solution [START_REF] Alamir | Stabilization of Nonlinear Systems Using Receding-Horizon Control Schemes: a Parametrized Approach for Fast Systems[END_REF], [START_REF] Van Donkelaar | Model predictive control with generalized input parametrization[END_REF], [START_REF] Muehlebach | Parametrized infinite-horizon model predictive control for linear time-invariant systems with input and state constraints[END_REF], [START_REF] Limaverde Filho | Trajectory tracking for a quadrotor system: A flatness-based nonlinear predictive control approach[END_REF]. The objective is to drastically reduce the number of control variables in the optimization problem while trying to keep the performance loss as small as possible. One of the major advantages of using a parameterization is the ability to decouple the prediction horizon from the control one. This is a very important improvement since for many systems the use of large prediction horizon is mandatory to ensure system stability in closed loop [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF]. However, in the classical approach, the two horizons are closely linked and increasing the prediction one necessarily implies enlarging the number of control variables. The complexity of the optimization problem, therefore, increases considerably. Even disregarding the issue of system instabilities, it is often necessary to select a large prediction horizon so as to obtain good closed-loop behavior. This is especially true when dealing with small sampling periods.

In this paper, we investigate the use of different input parameterizations in MPC to control non-linear systems. In particular, we review few existing input parameterizations coupled with a simple single-shooting scheme that can solve the MPC optimization using Sequential Quadratic Programming (SQP) quickly and while enforcing control input constraints such as bounds on the maximal command and its rate. Thanks to the parameterization, the sub-problems arising at each SQP iteration are low-dimensional and can thus be solved extremely quickly. Indeed, we show that using the parameterized approach it is possible to attain sub-millisecond performances despite a very dense timesampling and outperform existing state-of-the-art algorithms both on a laptop and on a Raspberry Pi device.

The remainder of this paper is organized as follows. In the next section, we firstly detail the adopted nonlinear MPC formulation and continue by recalling the studied control parameterizations. We propose a simple solver for the parameterized MPC optimization problem in section II-C and benchmark its performances in section III. In particular, we analyze how they change in relation to the chosen number of parameters and compare against state of the art solvers from an existing framework. Finally, the last Section reports our conclusions and proposes future perspectives.

II. NONLINEAR MODEL PREDICTIVE CONTROL

In this section, we recall the formulation of Direct Shooting methods for constrained Nonlinear MPC, which is at the basis of the employed solver. We then discuss how the introduction of a parameterization can help in reducing the complexity of the described problem and recall some existing parameterizations. We finally propose an algorithm to compute optimal parameters, which is a SQP strategy based on the Gauss-Newton Hessian approximation.

A. Direct Shooting Nonlinear MPC

The problem that Model Predictive Control tries to solve can be stated informally as follows: given an evolution model of the controlled system, find an optimal control trajectory which can steer the current state of a system towards a desired configuration, while meeting feasibility requirements such as actuation limits.

A variety of mathematical formulations can be found in the literature [START_REF] Binder | Introduction to model based optimization of chemical processes on moving horizons[END_REF], [START_REF] Diehl | Fast direct multiple shooting algorithms for optimal robot control[END_REF], which can be classified in three main categories: (1) indirect methods, (2) dynamic programming and (3) direct approaches. Indirect methods formulate the control problem as a functional optimization that is solved using variational calculus, while dynamic programming formulations exploit Bellman's principle of optimality to obtain a control sequence via a recursive procedure. The solutions to these problems are hard to find in the general case, or can become computationally expensive due to the so-called "curse of dimensionality" [START_REF] Findeisen | An introduction to nonlinear model predictive control[END_REF].

Direct methods work instead by transcribing the problem as a finite-dimensional nonlinear optimization. First of all, the controlled system is modeled as a discrete-time one, characterized by a transition function f that, given the state x k ∈ R n and control u k ∈ R m at the discrete time step k, returns the corresponding state that should be reached by the system at the next time step k + 1, i.e.,

x k+1 = f (x k , u k ) (1)
The "informal" optimal control problem stated in the beginning of this section is then commonly written as the following nonlinear optimization:

min Np k=1 x k -x k 2 Q k + Np-1 k=0 u k -u k 2 R k (2a)
subject to:

x k+1 -f (x k , u k ) = 0 ∀k = 0, • • • , N p -1 (2b) c x (x k ) ≤ 0 ∀k = 1, • • • , N p (2c) c u (u k ) ≤ 0 ∀k = 0, • • • , N p -1 (2d)
wherein N p is known as prediction horizon, while x k and u k represent the desired state and control inputs at the discrete time step k. v A is used to denote the norm of vector v weighted by the positive semidefinite matrix A, such that v 2 A = v T Av. The constraints (2b) are introduced to ensure that the sequence of states x = x 1 , • • • , x Np is coherent with the sequence of control inputs u = u 0 , • • • , u Np-1 according to the discrete time evolution model of the system. Constraints (2c) and (2d) can optionally be added to enforce explicitly, by the functions c x and c u respectively, feasibility of the states and controls contained in x and u. Finally, it must be noted that in the problem above the value x 0 is assumed to be a known constant, corresponding to the actual state of the system when performing the optimization.

Model constraints (2b) can be enforced in two ways. The first one is to use both u and x as decision variables and to explicitly deal with these constraints during the optimization. This leads to a sub-family of direct methods known as multiple-shooting techniques [START_REF] Bock | A multiple shooting algorithm for direct solution of optimal control problems[END_REF], [START_REF] Albersmeyer | The lifted newton method and its application in optimization[END_REF]. The second option consists in using a single-shooting method, in which only the control input sequence u is used as free variables for the search. Predicted future states can then be evaluated iteratively starting from x 0 and the selected control inputs. Model constraints are thus implicitly dealt with by this evaluation process. In this paper, we have focused our attention on single-shooting strategies and therefore (2b) are assumed to be always satisfied during the optimization.

To obtain a control sequence given the current state x 0 , a fairly popular choice is to solve the optimization using Sequential Quadratic Programming with Gauss-Newton approximations [START_REF] Giftthaler | A family of iterative gauss-newton shooting methods for nonlinear optimal control[END_REF] of the Hessian of the objective in (2a). This iterative solution technique works by firstly approximating the original problem around the current guess of u with a linearly constrained quadratic minimization. Solving the quadratic sub-problem allows to obtain a variation δ u for the control input sequence, and a line search is then performed to form a new guess u for the original nonlinear optimization. The new solution takes the form u = u+sδ u , with s ∈ [0, 1] chosen so that it causes a sufficient decrease in a given merit function that penalizes constraints violation and favors decreases in the original objective [START_REF] Han | A globally convergent method for nonlinear programming[END_REF], [START_REF] Powell | A fast algorithm for nonlinearly constrained optimization calculations[END_REF]. The whole process (solution of the sub-problem and evaluation of a new candidate) is repeated until termination criteria are met -typically, until the step size is small enough or a maximum period of time has elapsed.

B. Complexity Reduction via Parameterization

In the optimization scheme presented above, performances highly depend on the prediction horizon. Ideally, one would wish to sample the system at a high frequency to reduce the inaccuracies introduced by discretization. At the same time, predicting over a longer period of time should lead to better results since it enables the algorithm to intelligently consider maneuvers that are globally optimal and steer the system to the desired state. It is thereby beneficial to increase the value of N p in both cases. However, the number of decision variables (and therefore the complexity of the nonlinear optimization) grows with the prediction horizon as well.

Since the optimization has to be performed at each iteration, it is fundamental to be able to converge to a solution quickly enough to meet real-time feasibility. This is particularly important when the controlled system is characterized by a fast dynamics, presents instability or features limited computational power. In order to reduce the complexity of the problem to be solved, different strategies can be employed that aims at decreasing the number of decision variables.

One common practice consists in allowing only the first N c ≤ N p control samples to be freely allocated in the optimization [START_REF] Allibert | Predictive control for constrained image-based visual servoing[END_REF], [START_REF] Perez-Morales | Multi-sensor-based predictive control for autonomous parking in presence of pedestrians[END_REF]. The remaining N p -N c samples are instead set to the same value of u Nc-1 , as illustrated in fig. 1(a). This allows the number of decision variables to be substantially reduced, while still being able to foresee a sufficiently large amount of future samples. Nonetheless, concentrating all the "degrees of freedom" of the optimization in the beginning of the prediction interval can lead to sub-optimal performances. Since the last value of the control sequence is repeated a large number of times, its influence on the final state reached by the system is intuitively much larger than that of the first samples. Furthermore, if the system presents instabilities, constant signals can lead to significant displacements. The control action in the beginning the prediction horizon might thus be required to preemptively compensate for such large displacement, possibly leading to unnecessarily large motions and less regular control actions.

For the reasons above, it would be desirable to distribute the degrees of freedom of the optimization along the prediction horizon. One possible way to do so could consist in selecting, as "free samples", values from u that are not all at the beginning of the prediction horizon, but rather uniformly spread along it. Remaining values could instead be obtained by "holding" the free samples constant until the next one is reached (as shown in fig. 1(b)) or alternatively using linear interpolation (fig. 1(c)). These strategies are known in the literature as Move Blocking [START_REF] Cagienard | Move blocking strategies in receding horizon control[END_REF].

A further, more general, option can be considered to reduce the number of optimization variables, which relies on the choice of a set of N b basis functions

φ i : R + → R (i = 1, • • • , N b ).
The control sequence can then be generated from these functions by linear combination as:

u k = N b i=1 φ i (k∆t)η i (3)
where the combination coefficients η i are to be determined by the optimization and ∆t is the sampling time used in the prediction model. Some examples in this sense are damped Laguerre polynomials [START_REF] Muehlebach | Parametrized infinite-horizon model predictive control for linear time-invariant systems with input and state constraints[END_REF] and Haar wavelets [START_REF] Lee | Improving computational efficiency of model predictive control algorithm using wavelet transformation[END_REF]. In this paper, we propose to use a set of exponentially damped polynomials in the form: All functions from this basis have a unique maximum in R + located at t = iτ , such that φ i (iτ ) = 1. It is thus reasonably easy to tune the shape constant τ , which simply controls the locations of the maxima. To the best of our knowledge, this is the first time such basis is employed in a control parameterization.

φ i (t) = et iτ i e -t/τ (4)
The interest in using a basis of functions is that the overall control sequence can "inherit" some desirable properties from its generating functions. As an example, if all φ i are smooth (such as in ( 4)) then the generated sequences naturally feature reduced variations between consecutive samples. In addition, changes in a single combination parameter usually affects a whole portion of the control sequence, rather than a specific one. From a practical point of view, this implies that even with a reduced number of variables it is possible to generate a rather wide variety of control profiles.

All the strategies described so far can be unified under the same concept of control parameterization, which is a mapping from a given vector of parameters η ∈ R Nη to the control sequence u. Furthermore, all the parameterizations discussed above can be shown to be linear in η, i.e.,

u =   u 0 • • • u Np-1   =   Π 0 • • • Π Np-1   η = Πη (5)
In particular, it can be shown that for the first three parameterizations, which will be referred to as "Simple", "ZOH" (Zero-Order-Holder) and "LERP" (linear interpolation), the parameters η correspond to the values of the free control samples (see fig. 1) and that their combination matrices Π have a special banded form. On the contrary, function basis parameterizations feature dense combination matrices:

Π =      φ 1 (0) • • • φ N b (0) φ 1 (∆t) • • • φ N b (∆t) . . . . . . . . . φ 1 (∆t(N p -1)) • • • φ N b (∆t(N p -1))      (6)
with the parameters η being the stack of vectors η i from (3).

C. SQP Solver for Parameterized MPC

In this section we detail a Sequential Quadratic Programming solver that can be employed for the solution of the MPC optimization problem (2) when using parameterizations. The proposed algorithm works by formulating a quadratic subproblem, whose solution allows to update the current guess of the parameters η. These two steps (solution of the subproblem and update of the parameters) is repeated until termination conditions are met. As mentioned before, we consider in this paper a single-shooting strategy in which (2b) is directly satisfied, while concerning other constraints we pose limits on the magnitude of the control actions and its variation, i.e.,

u min ≤ u k ≤ u max ∀k = 0, • • • , N p -1 (7a) -∆u ≤ u k -u k-1 ≤ ∆u ∀k = 0, • • • , N p -1 (7b)
Concerning the definition of the quadratic sub-problem, we firstly consider a current guess η of the parameters, from which the control sequence u is computed using [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF]. Afterwards, predicted states x are evaluated recursively as function of the initial state x 0 and u using the discrete-time model of the system. States and controls in the objective are then linearized around η + δ η , with δ η representing variation of the current parameters. Injecting the linearized states into the original objective (2a) leads to

min 1 2 δ T η Hδ η + g T δ η (8a)
with

H = Np k=1 ∂x k ∂η T Q k ∂x k ∂η + Np-1 k=0 Π T k R k Π k (8b) g T = Np k=1 (x k -x k ) T Q k ∂x k ∂η + Np-1 k=0 (u k -u k ) T R k Π k (8c)
Note that H corresponds to a Newton-approximation of the Hessian of the original objective function. In the relations above, the Jacobians of the states with respect to the parameters can be evaluated using the recursive relation

∂x k+1 ∂η = ∂f ∂x ∂x k ∂η + ∂f ∂u Π k (9)
wherein the derivatives of f are computed at the k-th state and control samples.

Constraints [START_REF] Van Donkelaar | Model predictive control with generalized input parametrization[END_REF] can be easily translated into constraints on the parameters variation δ η thanks to [START_REF] Bemporad | Model predictive control based on linear programming -the explicit solution[END_REF]. However, insights on the used parameterizations can lead to more efficient solutions. First of all, in simple, ZOH and LERP parameterizations the parameters represent a set of control samples, and it is easy to understand that other values cannot exceed (be less than) the largest (smallest) free samples. As a consequence, control bounds can be readily translated into parameter limits:

η min ≤ η + δ η ≤ η max (10)
Similarly, variation constraints need not to be enforced at each discrete sample k, but rather at the free samples j only, giving:

-∆u ≤ (Π j -Π j-1 ) (η + δ η ) ≤ ∆u (11)
When function bases are used to parameterize the control inputs, in the general case there is no way to reduce the number of constraints [START_REF] Van Donkelaar | Model predictive control with generalized input parametrization[END_REF]. However, if the chosen functions are smooth, it is reasonable to assume that values of consecutive samples will not be too dissimilar. In particular, if a control sample located at a given discrete instant j is within bounds, chances are that neighboring samples will either be within bounds or exceed control limits by a small margin. Furthermore, as the inputs are generated by differentiable functions, control variations u k -u k-1 can be approximated as uk ∆t. Once again, under smoothness assumption, it is reasonable to reduce the number of points at which variation bounds are imposed. To summarize, the set of valid parameters can be approximated by

u min -∆u ≤ Π j ∆tΠ d j (η + δ η ) ≤ u max ∆u (12a)
which are applied only at some instants j selected by the user. In the relations above, matrix Π d j is defined in terms of the derivatives of the basis functions as:

Π d j = φ1 (j∆t) • • • φN b (j∆t) (13)
The sub-problem defined by (8a) and the constraints ( 7) is a dense QP optimization that can be solved by state-of-theart solvers, such as qpOASES [START_REF] Ferreau | qpoases: A parametric active-set algorithm for quadratic programming[END_REF], HPIPM [START_REF] Frison | HPIPM: a high-performance quadratic programming framework for model predictive control[END_REF] or OSQP [START_REF] Stellato | OSQP: an operator splitting solver for quadratic programs[END_REF]. Once the solution is computed, a line search algorithm can be used to update the current guess of the parameters as η = η + sδ η with s ∈ [0, 1]. The algorithm can then continue by repeating the previous steps over and over until convergence.

It must be noted that, in a classical direct single-shooting formulation, the number of decision variables, control bounds and control variation constraints are all proportional to the prediction horizon N p . Instead, in the parameterized approach proposed here, they are all proportional to N η . As shown in the next section, the use of a reduced number of parameters allows for considerable speed-ups in the solution process, without degrading the control performance.

III. SIMULATIONS

In this section we benchmark an implementation of the proposed solver in a simulated environment. The parameterization and algorithms detailed in previous sections were implemented as a C++ library 1 , using Eigen [START_REF] Guennebaud | Eigen v3[END_REF] for linear algebra and qpOASES [START_REF] Ferreau | qpoases: A parametric active-set algorithm for quadratic programming[END_REF] to solve the dense QP sub-problems. Different parameterizations with increasing number of parameters have been tested to determine the performances of the solver and a comparison with state of the art methods from the acados framework [START_REF] Verschueren | acados-a modular open-source framework for fast embedded optimal control[END_REF] is included.

The system to be controlled is an inverted pendulum with moving base (sometimes referred to as cart-pole system), which is often used to investigate new approaches in nonlinear control theory. The objective here is not to propose a new control strategy for this type of machine but simply to compare the efficiency of the proposed approach with the state of the art in MPC on a simple system. The state is given by the quadruplet (p, θ, v, ω), representing respectively the position of the cart, the angle of the pendulum (θ = 0 meaning that the pendulum points downward), the linear velocity of the base and the angular rate of the pendulum. The control input of the system is a horizontal force pushing the cart. The continuous-time model of the system can be obtained using a Lagrangian approach, which gives the linear and angular accelerations as function of the state and input:

p = I 0 F + µ θ2 sin θ + gµ 2 sin θ cos θ m t I 0 -µ 2 cos 2 θ (14a) θ = -F + µ θ2 sin θ µ cos θ -m t gµ sin θ m t I 0 -µ 2 cos 2 θ (14b)
wherein g is the gravity acceleration constant, m t the total mass of the system (base and pendulum), I 0 the moment of inertia of the pendulum with respect to the pivot point and µ the product of the mass of the pendulum times the distance between its center of mass and the pivot. Numerical values used in the simulations are: g = 9.806 m s -2 , m t = 0.3678 kg, I 0 = 0.011 852 kg m 2 and µ = 0.035 34 kg m. Different sets of simulations are presented in the sequel: first, we focus on a task in which the pendulum has to be kept in balance while also tracking a time-varying trajectory with the base. In this first set of simulations, all parameterizations presented in previous sections are tested multiple times with different values of N η . Afterwards, the results obtained with two parameterizations with N η = 5 are compared to four algorithms from acados. Finally, to test the effectiveness of the parameterized approach in a more challenging scenario, a swing-up task is considered.

In all simulations, the time horizon used for the predictions was set to 1 s, and a very high control frequency was employed to test the algorithms under challenging conditions. In particular, all presented simulations were performed with both ∆t = 10 ms and ∆t = 2 ms, leading to a prediction horizon of N p = 100 and N p = 500 samples respectively. Tests were run on two different platforms: (1) a Laptop 1 We are polishing the code and releasing it on GitHub in the final version.

computer featuring an Intel Core i7-8750H @2.2 GHz CPU and 16GB of RAM; (2) a Raspberry Pi 400 featuring a Broadcom BCM2711 Cortex-A72 (ARM v8) @1.8 GHz processor and 4GB of RAM. In this way, it was possible to get some insights on the performances of the implemented solver both on a workstation and an embedded device.

A. Performance benchmarking

In this first set of tests, the objective was to gather insight about the performances of parameterizations for an increasing number of parameters. Each simulation lasts 30 s, with the pendulum starting in the state (p, θ, v, ω) = (0, π, 0, 0) and having to track the desired states (p (t), π, 0, 0), with the position defined as:

p (t) =      -0.5 t < 7.5 0.5 7.5 ≤ t < 15 0.2 sin (0.4πt) t ≥ 15 (15)
Weighting matrices in [START_REF] Muehlebach | Parametrized infinite-horizon model predictive control for linear time-invariant systems with input and state constraints[END_REF] were defined as Q k = diag(200, 50, 7, 2)/N p and R k = 1/N p , and only control bounds (7a) were included, with u max = -u min = 0.1 N.

Concerning the tested parameterizations, the simple, ZOH and LERP ones were tested for N η = 2, • • • , 20, with free samples spread uniformly along the prediction horizon in the case of ZOH and LERP parameterizations. In addition, a further parameterization using a basis of exponentially damped polynomials was tested (later referred to as "Poly", with N b = N η = 2, • • • , 10 parameters. In each of these cases, the parameter τ was set to 1/(N η -2), so that the maxima of the basis functions are regularly spread along the prediction horizon. In addition, control bounds (7a) were imposed on 2N η samples, uniformly distributed along the prediction horizon.

For each parameterization, prediction horizon and number of parameters, the same simulation was repeated 10 times to be able to obtain statistical data on the amount of time required to solve the nonlinear optimization problems. Results are shown in figs. 2 and 3 for experiments run on a laptop and a Raspberry Pi respectively. For each parameterization, two types of graphs are included: average run-times (per iteration) in microseconds vs the number of parameters (figs. The results show that the simple parameterization consistently requires less computation time than other strategies, but that it fails to converge to a good solution, with tracking errors that remain more or less the same independently from the number of parameters used. This can be justified by the fact that parameters are able to affect only a small portion of the prediction horizon, with most of the control inputs being handled by a single parameter. ZOH and LERP parameterizations have very similar performances in terms of optimization times, with ZOH usually taking slightly more than LERP. However, the former features worse tracking performances for an equal number of parameters. Poly parameterization takes the longest computation times, which can be explained by the fact that computing the control sequence from parameters requires more flops and that the number of constraints is larger than in other cases. Nonetheless, it generally reaches good tracking performances with less computation times than the ZOH parameterization.

It is interesting to notice that in multiple cases the average optimization times are less than a millisecond. In particular, the results show that almost all the parameterizations could have been employed in real-time experiments in the case of a laptop. The average optimization times on the Raspberry board are larger, but still compatible with real-time requirements in the case of N p = 100.

B. Comparison against existing solvers: trajectory tracking

To better place the obtained results within the state of the art, the same simulations performed in the previous section were implemented also using the tools contained in acados. All tuning parameters of the problem, such as the weights in the objective function, were taken as in the previous section to be as fair as possible in the comparison. Termination conditions were also chosen to ensure that the results from the different strategies had the same level of optimality. Four multiple-shooting solvers were tested: one based on qpOASES and three HPIPM variants featuring full, partial and no condensing. As before, the simulations were repeated 10 times per solver, with N p set to 100 or 500.

Average time performances (with 95 % confidence intervals) at each iteration of the simulation are plotted in fig. 4 for all acados' algorithms and for LERP and Poly parameterizations2 (with N η = 5). In addition, for the case N p = 100, we include the results obtained using a simple parameterization with N η = N p = 100. This corresponds to a MPC problem with full degrees of freedom and is therefore labeled as "full" in the results. In addition, a of numerical values is included in the top-half of table I. In particular, for each algorithm we report the average optimization time over all iterations and all simulations in the different conditions. In addition, we evaluated a numerical index that quantifies how faster an algorithm is with respect to a given "baseline". This index is evaluated as:

speed gain = Niters i=1 T bl opt [i] T solver opt [i]
N iters [START_REF] Han | A globally convergent method for nonlinear programming[END_REF] where N iters is the number of iterations per simulation and T solver opt

[i] and T bl opt [i] are the average optimization times required at the i-th iteration to solve the SQP problem using respectively the given solver and the baseline (therefore, it does not necessarily coincide with the overall average time of the baseline divided by that of the solver).

Two observations are due. First of all, the two parameterized approaches consistently take shorter computation times than acados' solvers (they are always at the bottom of the plots shown in fig. 4). The Poly and in particular the LERP parameterizations both converge to a good solution faster than every solver from acados -especially on the Raspberry board, where they run more than 10 times faster on average. Furthermore, even though the values have been omitted for brevity, the average tracking errors obtained with these two parameterizations are slightly smaller than those of acados' algorithms, showing that shorter optimization times have not been obtained by sacrificing optimality. The second important remark is that the full parameterization takes much longer than any algorithm in acados. Since its working principles are close to those of the qpOASESbased solver in acados, this suggests that there is possibly room for further improvement in the efficiency of our software implementation.

C. Comparison against existing solvers: swing-up

The last set of simulations considered here involves a swing-up task, with the pendulum starting from the origin of its state space and having to reach (0, π, 0, 0) as desired configuration. Weighting matrices were changed to Q k = diag(10 3 , 10 3 , 10 -2 , 10 -2 )/N p and R k = 10 -2 /N p , and the control bounds relaxed to u max = -u min = 12 N. To make the problem more challenging, we also added variation constraints 3 (7b) with ∆u = 3∆tu max .

Considering the proposed solver, we decided to find for each parameterization the minimum number of N η that was sufficient to complete the task within 3 s. It turns out that the simple and ZOH parameterizations require at least 16 and 33 parameters respectively (see fig. 1). Instead, LERP and Poly parameterizations were successful with as few as 5 and 3 parameters each. Moreover, the full parameterization and the qpOASES-based solver in acados both failed to converge in a reasonable amount of time and were thus not included in the analysis.

Similarly to what done in the previous section, we show plots of the average optimization times during the simulations in fig. 5, while numerical results are included in the bottom-half of table I. Once again, the parameterized 3 Since in acados this is not directly possible, we added the force as a new state variable and used its derivative as the new control. In this way, both the actuation force (as a state) and its variation (being the new input) can be bounded. methods LERP and Poly are confirmed to be faster than existing solvers, by quite a large margin. In particular, it is worth remarking that on the Raspberry Pi board with N p = 100, only the parameterized approaches Simple, LERP and Poly converged within the sampling period of ∆t = 10 ms for more than half of the simulation.

IV. CONCLUSIONS

In this work we reviewed the benefits of using input parameterizations in Model Predictive Control to achieve fast and reliable performances. The advantages of these strategies are twofold. On one hand, they allow to reduce the number of optimization variables and of the constraints, making the problem solvable even on embedded hardware with very limited computational resources. On the other, they improve the performances of the control with respect to receding horizon techniques characterized by a small control horizon. We implemented a single-shooting SQP algorithm to solve the nonlinear MPC optimization and showed that, coupled with LERP and function basis parameterizations, it can consistently outperform state-of-the-art solvers.

The results obtained so far are very promising, and we believe that it is worth investigating parameterizations more in depth. In particular, we would like to focus our attention on different types of basis functions and extend the analysis to nonlinear parameterizations. The stability study of these parameterized problems is also an interesting track, that could be addressed by extending the stability results already obtained in the framework of the full parameterization [START_REF] Alamir | Stabilization of Nonlinear Systems Using Receding-Horizon Control Schemes: a Parametrized Approach for Fast Systems[END_REF].

Several works in the field suggest that multiple-shooting methods can provide more stable and fast solutions when formulating and solving the optimization instances. Indeed, ad-hoc condensing strategies for ZOH parameterizations have already been proposed [START_REF] Chen | Efficient move blocking strategy for multiple shooting-based non-linear model predictive control[END_REF]. We are confident that other parameterized approaches would also benefit from the use of lifted techniques. Furthermore, we would like to

Fig. 1 .

 1 Fig. 1. Example of control sequences u with a reduced number of degrees of freedom. Even if the prediction horizon is Np = 10, only 4 variables are freely assigned, while the remaining ones are defined in function of the free samples.

  2(a), 2(c), 3(a) and 3(c)) and average runtimes vs average position tracking errors (figs. 2(b), 2(d), 3(b) and 3(d)).

FinallyFig. 2 .

 2 Fig. 2. Performances (on a Laptop) of the proposed parameterizations with increasing number of parameters. In (a) and (c) the average optimization times are reported for increasing numbers of parameters, while (b) and (d) show the relation between optimization times and tracking errors (each curve corresponds to a parameterization, and each data-point to a specific number of parameters). The prediction horizon for (a) and (b) is Np = 100, while for (c) and (d) it is Np = 500.

Fig. 3 .

 3 Fig. 3. Performances of the proposed parameterizations on a Raspberry Pi board. The meaning of each plot is analogous to those depicted in fig. 2.

Fig. 4 .

 4 Fig. 4. Optimization times for the tracking task. HPIPM (f.c.), HPIPM (p.c), HPIPM: acados solvers based on HPIPM (with full, partial and no condensing); qpOASES: acados solver based on qpOASES; lerp-5, poly-5: proposed solver with Nη = 5. full: proposed solver with Nη = Np.

Fig. 5 .

 5 Fig. 5. Optimization times for the swing-up task.

Simple and ZOH were not considered here since with 5 parameters they did not converge to satisfactory solutions in terms of position error.

optimize our software library and add new features such as automatic differentiation and code generation, and possibly incorporating it into existing frameworks [START_REF] Giftthaler | The control toolbox -an open-source C++ library for robotics, optimal and model predictive control[END_REF], [START_REF] Verschueren | acados-a modular open-source framework for fast embedded optimal control[END_REF].


Egon Heuson

Augusto Etchegaray

Stephanie Léopoldine Filipe

Daniel Beretta

Mickaël Chevalier

Vincent Phalip

Dr François Coutte
email: francois.coutte@polytech-lille.fr

Screening of lipopeptide producing strains of Bacillus sp. using a new automated and sensitive fluorescence detection method †

Keywords: high-throughput screening, lipopeptides, surfactin, fluorescence, Bacillus CPC, cetylpyridinium chloride, FL, fluorescein, BTB, bromothymol blue, RP-UPLC, reversed phase ultra-performance liquid chromatography, MALDI-ToF, Matrix Assisted Laser Desorption Ionisation -Time of Flight. MIC: Minimum Inhibitory Concentration

Lipopeptides, such as surfactins are important biosurfactants produced by Bacillus sp. that find applications in many areas (environment, medicine and food industries).

Giving their importance, the use of simple detection methods will facilitate screening and quantification. In the present work, we described a completely automated workflow for the screening of lipopeptide producing strains and its quantification. Firstly, isolated colonies from environmental samples were automatically picked and inoculated in 96 wells growth plate. After overnight incubation, surfactin produced in the broth was quantified, using a new sensitive fluorescent method. The method uses fluorescein (FL), which is an anionic dye at neutral to alkaline pH and forms a stable complex with the cationic surfactant cetylpiridinium chloride (CPC), quenching fluorescence. Upon addition of surfactin or other lipopeptides, fluorescein is released from the CPC-FL complex and quantified. The robustness of this method was assessed by comparing the quantification results to those conventionally measured by RP-UPLC and the results of strain screening were confirmed by MALDI-ToF analysis. We report for the first time the successful application of this analytical method for high-throughput screening of novel lipopeptideproducing strains.

Highlights

Previous work based on colorimetric detection of surfactin uses the dye-complex CPC-BTB. In the present work, bromothymol blue was replaced by fluorescein, forming a stable complex at low concentration. By using automation, it was possible to design a complete workflow that is useful for quantification and high-throughput screening of novel surfactin-producing strains.

Biosurfactants are important secondary metabolites produced by certain bacterial strains. The importance of these compounds are their applications in environmental remediation [1], their antimicrobial activities [2,3] and synergistic action with other secondary metabolites [4]. One of the most studied and active biosurfactants are surfactins, an heptapeptide with strong surface activity [5,[START_REF] Sen | Biosurfactants[END_REF]. Two other compounds of the family of surfactins, namely pumilacidin and lichenysin have a structure very similar to surfactin. Other important lipopeptides produced by Bacillus sp. are iturins (iturin A, mycosubtilin and bacillomycin) and fengycins (fengycin A and B, plipastatin), which have antifungal activities. The structural features of surfactins include the presence of a betahydroxy fatty acid (C 12 -C 16 ) and the aminoacids L-glutamate and L-aspartate, which contribute with negative charges for this anionic biosurfactant [START_REF] Jacques | [END_REF]8]. The interest of the scientific community for these biosurfactant molecules, such as lipopeptides and particularly for surfactin has been increasing steadily since the last 10 years. For example, the number of scientific articles referring to surfactins over this period is 975 and 3450 more generally on lipopeptides, i.e. almost one scientific article per day (query on Scopus database the October 23 rd , 2018). The growing interest in these molecules leads to an increased work of the scientific community for the optimization of their production either by genetic and metabolic engineering [8,9] or process engineering [10] but the limits of these approaches are being achieved. In addition, the society is increasingly in demand for natural products that do not originate from genetically modified organisms. In this context, search for new strain producing these molecules with interesting production capacities or which could be naturally optimized becomes an important issue.

Quantification of surfactants is difficult to carry out by high-throughput methods [START_REF] Biniarz | Screening Concepts, Characterization and Structural Analysis of Microbial-Derived Bioactive Lipopeptides: A Review[END_REF]. The work of [START_REF] -Y. Chen | [END_REF] describes however, a method that is based on the effect of surfactant on the solution's meniscus and its consequence on the visualization of a grid that is placed under a microtiter plate [START_REF] -Y. Chen | [END_REF].

The quantification of surfactin (or other lipopeptides from Bacillus sp.) usually is carried out by liquid chromatography (HPLC). However, it requires a clean-up step by solid phase extraction [8]. Other methods such as surface tension measurement [13], oil spreading [14], blood hemolysis [15] or drop collapsing [16] test have often been used [13]. Recently Yang et al. [17] presented a colorimetric procedure. The method uses the anionic dye bromotimol blue (BTB) as indicator and the cationic surfactant CPC as mediator. However, the dye complex is not stable. It only last a few hours before precipitating and losing color.

In addition, the limit of detection seems to be around 0.1 g/L [17].

Here, we present a complete automated workflow for screening of lipopeptides (such as surfactin) producing strains (Fig. 1). This workflow is divided in three steps: a first step for automated colony picking from agar medium, a second step for incubating the strains on microtiter plates containing production medium and a third step for quantifying the lipopeptides produced using fluorescent methods, thus inspired on the work described by Yang et al. [17]. The method developed and optimized here uses the anionic fluorescent dye FL as an indicator dye that was mixed with CPC as mediator, resulting in a more stable complex, with increased sensitivity. In addition, it produces better curves, as there is a gain in linearity, besides producing an automated method that could be used for highthroughput screening of lipopeptide production and especially for surfactins. This workflow has been successfully tested to screen novel isolated strains of Bacillus sp. and to quantify the production of surfactins. The method was validated by ultra-performance liquid chromatography (UPLC) and a coupled MALDI-ToF analysis confirmed the species of the novel isolated strains and the type of lipopeptides produced.

Materials and Methods

CPC-FL method optimization

Three steps of optimization were used to design this new quantification method.

All experiments were performed in Greiner® black 96-flat-well polystyrene plates, incubated at 37 ± 1 °C and stirred for 30 s before measurement.

Initially, the excitation and emission wavelength of FL were optimized using a multimode microplate reader (SpectraMax i3, Tecan). The focalization height was also optimized and found to be ideal at 0.5 mm above the bottom of the well. FL is a fluorescent dye with absorption and emission in the visible range, respectively around 494 and 512 nm [18]. The results obtained here have shown that the fluorescein-displacing reaction can be optimally monitored using an excitation wavelength of λ ex = 482 ± 9 nm and a detection at λ em = 528 ± 15 nm. PMT gain was set to the highest value, with 6 flashes per read.

Secondly, relative fluorescence units (RFU) detected were correlated with FL concentration using a calibration range from 0 to 0.01 mM in 0.1 M sodium phosphate buffer pH 8.0. Results are presented in Fig. 2a. The optimal CPC/FL ratio and CPC concentration were determined using a 0.045 to 0.360 mM of CPC with 0.009 mM of FL, in presence of 0.7 to 45 mg/L of surfactin final concentration in the reaction mix. Results are presented in Fig. 2b. Experiments were carried out 3 times and with a technical triplicate. The mean values and standard deviation are presented.

Determination of the MIC of CPC

The minimum inhibitory concentration (MIC) of CPC alone and the mixture of CPC-FL were determined in microtiter plate by cultivating the strain Bacillus subtilis ATCC21332 [ATCC Collection] with different concentration of CPC. The wells of a 96 well plate, containing 100 µL of LB medium with different concentration of CPC (from 0 to 0.2 mM) were inoculated with 10 µL of pre-culture to obtain an initial optical density of 0.1.

After 24 hours of incubation at 37°C under 200 rpm orbital shaking, the growth in each well was measured by optical density at 605 nm (SpectraMax i3, Molecular Devices).

Experiment was carried out in triplicate.

Screening of environmental strains using the CPC-FL method

Environmental samples were collected from a farm soil in the Haut de France region (France) and then freezed at -20°C before analysis. The samples were first diluted in sterile distilled water and heated at 80°C for 20 minutes, in order to select only sporulating strains. The samples were then streaked on 250*250 mm YPG agar plates containing 100 µg/mL amphotericin B and incubated at 37°C for 24 hours. Then, 10 isolated colonies were randomly selected regarding their different morphology, and robotically transplanted using the colony picker (QPix 460, Molecular Devices) into a deep-well plate (Bioblock) of 2.2 mL well capacity (Simport, Beloeil, Canada). The wellscontained 1 mL of modified Landy's medium that was adjusted to pH 7.0 using phosphate buffer as previously described [9] or 1 mL of LB medium. Landy medium, a chemically defined medium, is one of the most used medium for lipopeptide production. FL in modified Landy's and LB medium, which were previously incubated at 37 °C for 24 h at 200 rpm orbital shaking. Fluorescence measurements were performed after 30 s of orbital shaking at 30 °C, using λ ex = 482 ± 9 nm as excitation wavelength and λ em = 528 ± 15 nm for detection as described previously. Experiments were carried out 3 times and with a technical triplicate. The mean values and standard deviation are presented.

Lipopeptides quantification using RP-UPLC

Lipopeptides (i.e. iturins, fengycins and surfactins) were quantified using an ultraperformance liquid chromatography (UPLC) ACQUITY H-Class system from WATERS Waters PDA, Empower Software) of each peak were used to identify and quantify the eluted molecules.

Strain identification by MALDI-TOF-MS

Isolated bacterial strains from LB agar plate were loaded three times onto ground steel MALDI target according to the manufacturer's instructions (Bruker Daltonics, Bremen, Germany). Briefly, 3 bacterial colonies were taken off from agar plates with a 10 μL pipet tip, smeared onto the target and dried at room temperature. The on-target deposits were overlaid with 1 µL of 70% formic acid solution, dried at room temperature, overlaid again with 1 µL of matrix solution [10 mg/mL of HCCA dissolved in ACN/water/TFA (50/47.5/2.5; v/v/v)] and dried again [19,20]. Samples were automatically analyzed using MALDI-TOF mass-spectrometer (AutoflexSpeed from Bruker Daltonics) running Flexcontrol 3.4 software as described recently [19].

Results and discussion

In this work, we designed, optimized and assessed a fully automated workflow for screening bacterial lipopeptide (especially surfactin) production and quantification on microtiter plates using CPC-FL. The method is based on the release of FL upon interaction with anionic lipopeptides.

FL is a very useful dye that finds many applications in analytical chemistry [21].

The spectroscopic properties of FL are based on its high fluorescence emission at 520 nm for the deprotonated species [22]. In this work, we have explored the properties given by the complex formation between FL and CPC, which quenches fluorescence. The quenching phenomenon is potentially given by electrostatic interactions between the cationic surfactant and deprotonated FL at pH 8.0. At 25°C, the critical micellar concentration (CMC) of CPC in water is approximately 1 mM [23]. The working concentration of CPC used in the experiment was (0.09 mM), which is approximately 10 times below the CMC of CPC. Therefore, we expect that the interaction between CPC and FL is given by charge interactions occurring at the molecular level. However, the ion-pair formed between CPC and FL may decrease the solubility of CPC, by charge neutralization. Thus, it may induce surfactant self-aggregation. A similar pattern was observed for the interaction between CPC and another anionic dye (phenol red) [24]. For surfactin, the lowest concentrations of the analytical curve are close to its reported CMC (0.0075 mM) [25]. Therefore, it is possible that the interaction of surfactin with CPC involves adsorption of CPC on surfactin micelles. The binding of CPC to surfactin micelles is certainly stronger when compared to its binding to FL. The CPC-surfactin complex involves both electrostatic and hydrophobic interactions, possibly forming mixed micelles. Also, it must be considered that surfactin has two negative charges in the molecule, while fluorescein has only one negative charge at pH 8.0 [22]. Because of surfactin binding, free FL is released, and fluorescence is detected.

Several calibration curves were prepared to find the best parameters for the CPC-FL method to evaluate surfactin production, during the growth of Bacillus sp. Relative fluorescence units (RFU), which were detected in correlation with FL concentration, show good linearity, up to 0.01 mM (r2 = 0.9976) (Fig. 2a) and an excellent LOD (Limit Of Detection) of 1.5 x 10-4 mM, which was determined according to IUPAC specifications.

Considering linearity as an important parameter, the optimal concentration of CPC required to obtain the calibration curve was calculated. A 10-fold factor between CPC and FL concentration was found ideal, with a final concentration of 0.09 mM in CPC and 0.009 mM in FL, as shown on Fig. 2b. It is important to note that the parent method using BTB-CPC as the reporting dye complex has a problem of stability. In our experiments, the dye complex tends to precipitate. Spectrophotometric analysis showed a reduction of absorbance of 50% in 40 hours due to precipitation, while the FL-CPC dye remained stable for longer period, especially at the established concentration above (data not shown).

From this, the optimized method consists in the use of 0.009 mM FL, as indicator dye, mixed with 0.09 mM CPC in 0.09 M sodium phosphate buffer pH 8.0.

Using the above described conditions, the lipopeptides calibration curves in water (Fig. 2c), in modified Landy's medium (Fig. 2d) and in LB medium (Fig. 2e) were determined. A very low LOD (0.015 g/L) could be calculated for the calibration curve in water. The LOD for the modified Landy´s medium curve showed even better results (0.008 g/L) for surfactin, according to IUPAC specifications. In both cases, a perfect linearity was found up to 450 mg/L of surfactin (r 2 = 0.9912 and 0.9988 respectively). In LB, however, the LOD was found higher, and the linearity not as good as in the previous conditions, with respectively 0.044 g/L and r 2 = 0.9838. This is correlated to a much lower slope for the calibration curve in LB than in water and Landy, that might come from the very high quantities of peptides present in LB and that may interfere with the test. These samples were also analyzed by RP-UPLC and similar quantification results were obtained (data not shown). Better linearity was obtained using this CPC-FL method in comparison to CPC-BTB method (r 2 = 0.9793) [17]. For fengycin, the LOD found in water, Landy and LB media were of 0.041 g/L, 0.058 g/L and 0.051 g/L, respectively. Linearity coefficient were also acceptable with respectively r 2 = 0.9913, r 2 = 0.9985 and r 2 = 0.9747 in the three media.

Here again, the linearity was lower in the case of LB, showing the negative impact of this broth on the assay response. Finally, for mycosubtilin, LOD were of 0.096 g/L and 0.110 g/L for water and Landy respectively, with correlation coefficients of r 2 = 0.9967 and r 2 = 0.9952. No LOD could be calculated from the mycosubtilin calibration curve in LB as no linear correlation was observed between the concentrations. Interaction of iturinic compound with CPC is significantly less good because they are uncharged. Nevertheless, in a chemically defined medium, it can certainly interact via a hydrophobic interaction and can be detected. More globally, this shows again the good correlation existing between the slope and the LOD value.

The method was then evaluated in the complete workflow described above and schematically shown in Fig. 1 for the screening of unknow environmental sample and using 2 different reference strains of B. subtilis. The time required to pick and analyze the production capacity of, for example, 96 strains (which corresponds to a 96-well plate) using this automated workflow, was 2 h (compared to at least 35 min of UPLC analysis per sample). To this period, 16 h of incubation is required on Petri dishes for bacterial growth, in addition to 24 h for surfactin production, after inoculation of single colonies on the liquid broth. However, these two incubation periods are inherent to all methods. At the end of the workflow, the new CPC-FL quantification method was compared to the usual method using RP-UPLC. The results are presented in Fig. 3. During the first step of screening 10 different colonies were selected and picked but only 5 of these grew on LB and modified Landy's medium. Analysis of culture supernatants by RP-UPLC indicated no fengycins or iturins production by any of the analyzed strains, including the producing strain ATCC 21332. The low cultivation time (24h) has not been suitable for the production of these lipopeptides, which are most often produced during the stationary phase of growth [26]. RP-UPLC analysis of the 168 strain's fermentation broths (a nonproducer of lipopeptides) reveals no surfactin production while CPC-FL method gives a response of 0.0098 ± 0.005 g/L in LB medium and 0.039 ± 0.011 g/L in modified Landy's medium. This result cannot come from the medium since the non-inoculated medium (negative control) used in this experiment produces almost no response. Nevertheless, these background noise values were used to normalize (by subtraction) the results obtained for the other strains screened.

The production results obtained with the reference strain ATCC 21332 show the effectiveness of the detection method described herein. It will be noted a good match of quantifications in the LB medium at around 0.015 g/L while a larger difference is observed in the Landy's medium between the two methods of quantification. Indeed surfactin production was measured at 0.38 g/L by RP-UPLC and at 0.44 g/L by the CPC-FL method.

For the other five strains analyzed here, it can observe initially that two of these do not produce lipopeptides (strains A and E). For these two strains, the two quantification methods gave similar results, which are very close to the background noise obtained for the strain 168. This result is confirmed by the MALDI-ToF analysis, where no characteristic peak of the three families of lipopeptides (surfactins, iturins or fengycins) was detected (data not shown). For the other three strains, isolated from environmental sample, the production of surfactin was successfully quantified by both the CPC-FL and RP-UPLC to be in the range of 0.029 g/L to 0.053 g/L in the LB medium. The difference between the two methods of quantifications was not judged as statistically significant.

These results were also confirmed for strain B in the modified Landy's medium.

Nevertheless, for strains C and D, the quantification obtained by RP-UPLC after cultivation in the modified Landy's medium is very different from that obtained by the CPC-FL method. These results seem to reveal the production of other compounds which can interfere in the assay. One of the explanations would be an interference caused by primary metabolites resulting from the degradation of glucose that is present in Landy medium and not in the LB, such as acetate and lactate. These two compounds are very often produced by B. subtilis under low oxygen transfer, as is the case of liquid growth on microtiter plates [8,27]. To verify this, measurements were made with different concentrations of acetate and lactate (between 0 and 0.5 g/L) using the CPC-FL complex.

However, no fluorescein release was observed in the presence of these compounds (data not shown). Therefore, it is possible to infer that the difference between the UPLC results and the fluorescence-based quantification of surfactin, using a glucose based medium is due to production of additional primary or secondary metabolites by specific strains during cultivation. Interestingly, Yang et al., seem to have detected a similar pattern, when developing their method with CPC-BTB and using glucose based medium. This conclusion is based on the authors' specification that the optimal range for quantification, using the CPC-BTB method, is between 0.1 and 0.5 g/L [17]. This information suggests that other metabolites may give a corresponding false response at a concentration of 0.1 g/L. In addition to that, if longer incubation periods are to be analyzed, one must consider the contribution given by fengycins, which are additional anionic lipopeptide produced by Bacillus strains. In our studies, we found that this other biosurfactant give less sensitive response using the CPC-FL method, possibly because surfactin has 2 negative charges in comparison with fengycin. Finally, to validate our screening method further analysis was carried out by MALDI-ToF-MS. The analysis was performed directly on isolated colonies for species classification and also to determine their lipopeptide production pattern. For A novel analytical method was described that is very efficient to screen new lipopeptides producing strains. The method allows the detection and quantification of lipopeptides, such as surfactins, from environmental isolates in a single run. This completely automated method can be used also to collect kinetic data on surfactin production, study novel medium compositions or select novel biosurfactant producers, without requiring an expensive and time-consuming procedure such as liquid chromatography. The CPC-FL method was based on a previously described procedure using the dye BTB and it was shown to be more sensitive and present better linear response, besides it was adapted to a completely automated system and validated by RP-UPLC analysis. A parallel analysis of the CPC-FL detection method and RP-UPLC revealed the potentialities of the described method for screening novel surfactin-producing strains.

This was also corroborated using MALDI-ToF, which selects bacterial strains by their lipopeptide producing capacity. Nevertheless, special attention must be taken if medium highly enriched in a carbon source is used, which enables the cell to grow well and thus produce other molecules that can interfere in the assay. The use of a fluorescent reporting dye also opens opportunities for the CPC-FL method to be applied directly during the growth in micro-fermentation devices such as the Biolector® (from m2p-labs), which has inbuilt fluorescence detector using similar excitation and emission filters. To evaluate this, further studies will have to be conducted because CPC is toxic to microbial cells. Our studies also have shown that the MIC of CPC is approximately 0.1 mM.

Finally, the calibration

  curves of the different lipopeptides (surfactin, mycosubtilin and fengycin) were set by addition of 10 µL of a 7 to 450 mg/L lipopeptide solution in water, to 90 µL of 0.1 mM CPC and 0.01 mM FL solution in 0.1 M sodium phosphate buffer pH 8.0. Results are presented in Fig. 2c. Lipopeptides used for method optimization were supplied from Lipofabrik SAS (Villeneuve d'Ascq, France). The anionic fluorescent dye fluorescein (FL) (Uranin AP) and cetylpyridinium chloride (CPC) were purchased from Sigma Aldrich (Saint Louis, USA).

  Cultures were then set at 37°C during 24 hours under 200 rpm orbital shaking. Two reference strains of B. subtilis were also used, respectively as negative control, strain 168 a non-lipopeptides producer [Bacillus Stock Center Collection] [15] and positive control strain ATCC21332 a natural producer of surfactin and fengycin [ATCC collection]. The growth was first measured by optical density at 605 nm (SpectraMax i3, Molecular Devices) and strains were reisolated on separate LB agar plates for further analysis by MALDI-ToF. Then, 96 well-plates were centrifuged (Alegra X-30R, Beckman Coulter) for 30 min at 4 °C and 2700 g. A sample of 200 µL was kept at 4°C before RP-UPLC analysis. A sample of 10µL was transferred to a Greiner® black 96-flat-well plate. Sample preparation steps were performed on a Biomek FXp liquid handling station from Beckman Coulter. The fluorescent lipopeptide detection was automated on a multimode microplate reader (SpectraMax i3, Molecular Devices). The reaction mixture was composed of 10 μL of sample solution, and 90 μL of 0.1 mM CPC and 0.01 mM FL in 0.1 M sodium phosphate buffer pH 8.0. Calibration curves of each lipopeptide were prepared using 0.7 to 45 mg/L

(

  Saint-Quentin En Yvelines, France) and 2.1X50 mm BEH C18 Column. Lipopeptides were analysed using gradient with the initial solvent acetonitrile/water/trifluoroacetic acid 30:70:0.1 (vol/vol/vol) and reach 45:55:0.1 (vol/vol/vol) in 15 minutes to elute iturins. Then, a linear gradient to reach 55:45:0.1 (vol/vol/vol) in 10 min was applied. Under these conditions, fengycins are eluted. Following this step, another gradient of 5 minutes was performed to reach 80:20:0.1 (vol/vol/vol). Then surfactins were eluted isocratically during 5 min. Pure iturins A, fengycins and surfactins from Sigma-Aldrich were used as references. Retention time and second derivatives of UV-visible spectrum (Diode Array

2 . 6

 26 Mass spectra were processed using the BioTyper software (version 3.0; Bruker Daltonics) running with the BioTyper database version DB 5989, containing 5989 reference MALDI-TOF MS profiles (5298 of bacteria, 626 of yeasts and 65 of filamentous fungi). Matching between experimental MALDI-TOF MS profiles obtained from bacteria isolates and the reference MALDI-TOF MS profiles is expressed by BioTyper according to a Log (Score) and an associated colour code (green, yellow and red) as decribed recently [19]. Characterization of the lipopeptides production in novel isolates strains by MALDI-TOF-MS On the same deposits (smears) produced previously for strain identification, the MS signals were acquired in reflectron mode in the 700-3,500 Da m/z range by summing 8000 laser shot spectra, according to the manufacturer's automatic method RP_700-3500.par (voltage values of ion sources #1 and #2 set as 19 and 16.70 keV, respectively; voltage values of reflectron #1 and #2 set as 21 and 9.5 keV, respectively; lens tension = 6.9 keV; pulsed extraction = 120 nsec; laser intensity between 30 and 50%). Mass spectra were visualized using FlexAnalysis software (version 3.4; Bruker Daltonics).

  strains B, C and D, the characteristic peaks (m/z) of surfactin ([M + Na] + 1030, 1044, 1058, [M + K] + 1074) were found as shown for example for the strain B (Fig 4). Identification results also show that these strains belong to the species of Bacillus subtilis, Brevibacillus borstelensis, Bacillus licheniformis and Bacillus oleronius (data not shown). All these species belong to the Bacillus genus and are known to produce lipopeptides of the surfactins family such as surfactins, pumilacidins or lichenysins.

  Figure 1. Complete workflow for high-throughput screening of strains producing anionic lipopeptides such as surfactin. The workflow is divided in three steps: firstly colony picking; followed by incubation on microtiter plates and quantification based on fluorescein displacement from the CPC-FL complex.

Figure 2 .

 2 Figure 2. a) Fluorescein calibration curve. b) Relative fluorescence units measured in function of surfactin and CPC (✕ 0.045 mM; + 0.068 mM; ▲ 0.090 mM; ■ 0.135 mM; ◆

Figure 3 .

 3 Figure 3. Quantification of surfactin produced by the different strains after 24 h of growth in (a) LB medium and (b) modified Landy's medium at 37°C and under agitation (200 rpm) using CPC-FL method (black bars) and compared to RP-UPLC method (white bars). Background noise values obtained in the case of the 168 strain's fermentation broth (0.0098 ± 0.005 g/L in LB medium and 0.039 ± 0.011 g/L in modified Landy's medium) were used to normalize (by subtraction) the results obtained for the other strains. Mean values and standard deviation are presented (n=9).

Figure 4 .

 4 Figure 4. MALDI-ToF-MS analysis of the lipopeptides produced by the environmental strain B. The MS signals were acquired in reflectron mode in the 700-3,500 Da m/z range by summing 8000 laser shot spectra

Acknowledgements

The authors thank Anika Mros from School Centrale Lille for her technical assistance and the ALIBIOTECH program funding administered by the Hauts-de-France Region and the RealCat platform. The REALCAT platform is benefiting from a state subsidy administrated by the French National Research Agency (ANR) within the frame of the 'Future Investments' program (PIA), with the contractual reference 'ANR-11-EQPX-0037'. The European Union, through the ERDF funding administered by the Hauts-de-France Region, has co-financed the platform. Centrale Lille, the CNRS, and Lille 1 University as well as the Centrale Initiatives Foundation, are thanked for their financial contributions to the acquisition and implementation of the equipment of the REALCAT platform. The authors would also like to thank the Fundação de Amparo à Pesquisa de São Paulo (FAPESP), Brazil, respectively grants 13/20570-6 and 17/13330-0. Daniel Beretta was the recipient of a scholarship from PIBIC-CNPq, Brazil.

Conflict of Interest

The authors declare that they have no conflict of interest.