Robin Strudel

Alexander Pashevich

Igor Kalevatykh

Ivan Laptev

Josef Sivic

Cordelia Schmid

Learning to combine primitive skills: A step towards versatile robotic manipulation §

Manipulation tasks such as preparing a meal or assembling furniture remain highly challenging for robotics and vision. Traditional task and motion planning (TAMP) methods can solve complex tasks but require full state observability and are not adapted to dynamic scene changes. Recent learning methods can operate directly on visual inputs but typically require many demonstrations and/or task-specific reward engineering. In this work we aim to overcome previous limitations and propose a reinforcement learning (RL) approach to task planning that learns to combine primitive skills. First, compared to previous learning methods, our approach requires neither intermediate rewards nor complete task demonstrations during training. Second, we demonstrate the versatility of our vision-based task planning in challenging settings with temporary occlusions and dynamic scene changes. Third, we propose an efficient training of basic skills from few synthetic demonstrations by exploring recent CNN architectures and data augmentation. Notably, while all of our policies are learned on visual inputs in simulated environments, we demonstrate the successful transfer and high success rates when applying such policies to manipulation tasks on a real UR5 robotic arm.

I. INTRODUCTION

In this work we consider visually guided robotics manipulations and aim to learn robust visuomotor control policies for particular tasks. Autonomous manipulations such as assembling IKEA furniture [START_REF] Suarez-Ruiz | Can robots assemble an IKEA chair?[END_REF] remain highly challenging given the complexity of real environments as well as partial and uncertain observations provided by the sensors. Successful methods for task and motion planning (TAMP) [START_REF] Srivastava | Combined task and motion planning through an extensible plannerindependent interface layer[END_REF]- [START_REF] Toussaint | Logic-geometric programming: An optimization-based approach to combined task and motion planning[END_REF] achieve impressive results for complex tasks but often rely on limiting assumptions such as the full state observability and known 3D shape models for manipulated objects. Moreover, TAMP methods usually complete planning before execution and are not robust to dynamic scene changes.

Recent learning methods aim to learn visuomotor control policies directly from image inputs. Imitation learning (IL) [START_REF] Pomerleau | ALVINN: An autonomous land vehicle in a neural network[END_REF]- [START_REF] Ng | Algorithms for inverse reinforcement learning[END_REF] is a supervised approach that can be used to learn simple skills from expert demonstrations. One drawback of IL is its difficulty to handle new states that have not been observed during demonstrations. While increasing the number of demonstrations helps to alleviate this issue, an exhaustive sampling of action sequences and scenarios becomes impractical for long and complex tasks.

In contrast, reinforcement learning (RL) requires little supervision and achieves excellent results for some challenging tasks [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF], [START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF]. RL explores previously unseen scenarios and, hence, can generalize beyond expert demonstrations. As full exploration is exponentially hard and becomes impractical for problems with long horizons, RL often relies on careful engineering of rewards designed for specific tasks.

Common tasks such as preparing food or assembling furniture require long sequences of steps composed of many different actions. Such tasks have long horizons and, hence, are difficult to solve by either RL or IL methods alone. To address this issue, we propose a RL-based method that learns to combine simple imitation-based policies. Our approach simplifies RL by reducing its exploration to sequences with a limited number of primitive actions, that we call skills.

Given a set of pre-trained skills such as "grasp a cube" or "pour from a cup", we train RL with sparse binary rewards corresponding to the correct/incorrect execution of the full task. While hierarchical policies have been proposed in the past [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF], our approach can learn composite manipulations using no intermediate rewards and no demonstrations of full tasks. Hence, the proposed method can be directly applied to learn new tasks. See Figure 1 for an overview of our approach.

Our skills are low-level visuomotor controllers learned from synthetic demonstrated trajectories with behavioral cloning (BC) [START_REF] Pomerleau | ALVINN: An autonomous land vehicle in a neural network[END_REF]. Examples of skills include go to the bowl, grasp the object, pour from the held object, release the held object, etc. We automatically generate expert synthetic demonstrations and learn corresponding skills in simulated environments. We also minimize the number of required demonstrations by choosing appropriate CNN architectures and data augmentation methods. Our approach is shown to compare favorably to the state of the art [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF] on the FetchPickPlace test environment [START_REF] Plappert | Multi-goal reinforcement learning: Challenging robotics environments and request for research[END_REF]. Moreover, using recent techniques for domain adaptation [START_REF] Pashevich | Learning to augment synthetic images for sim2real policy transfer[END_REF] we demonstrate the successful transfer and high accuracy of our simulator-trained policies when tested on a real robot We compare our approach with two classical methods: (a) an open-loop controller estimating object positions and applying a standard motion planner (b) a closed-loop controller adapting the control to re-estimated object positions. We show the robustness of our approach to a variety of perturbations. The perturbations include dynamic change of object positions, new object instances and temporary object occlusions. The versatility of learned policies comes from both the reactivity of the BC learned skills and the ability of the RL master policy to re-plan in case of failure. Our approach allows to compute adaptive control and planning in real-time.

In summary, this work makes the following contributions. (i) We propose to learn robust RL policies that combine BC skills to solve composite tasks. (ii) We present sample efficient training of BC skills and demonstrate an improvement compared to the state of the art. (iii) We demonstrate successful learning of relatively complex manipulation tasks with neither intermediate rewards nor full demonstrations. (iv) We successfully transfer and execute policies learned in simulation to real robot setups. (v) We show successful task completion in the presence of perturbations.

Our simulation environments and the code used in this work is publicly available on our website [START_REF]Learning to combine primitive skills: A step towards versatile robotic manipulation[END_REF].

II. RELATED WORK

Our work is related to robotics manipulation such as grasping [START_REF] Lampe | Acquiring visual servoing reaching and grasping skills using neural reinforcement learning[END_REF], opening doors [17], screwing the cap of a bottle [START_REF] Levine | End-to-end training of deep visuomotor policies[END_REF] and cube stacking [START_REF] Popov | Data-efficient Deep Reinforcement Learning for Dexterous Manipulation[END_REF]. Such tasks have been addressed by various methods including imitation learning (IL) [START_REF] Duan | One-shot imitation learning[END_REF] and reinforcement learning (RL) [START_REF] Riedmiller | Learning by playing -Solving sparse reward tasks from scratch[END_REF].

Imitation learning (IL).

A neural network is trained to solve a task by observing demonstrations. Approaches include behavioral cloning (BC) [START_REF] Pomerleau | ALVINN: An autonomous land vehicle in a neural network[END_REF] and inverse reinforcement learning [START_REF] Ng | Algorithms for inverse reinforcement learning[END_REF]. BC learns a function that maps states to expert actions [START_REF] Ross | Reinforcement and imitation learning via interactive no-regret learning[END_REF], [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF], whereas inverse reinforcement learning learns a reward function from demonstrations in order to solve the task with RL [START_REF] Popov | Data-efficient Deep Reinforcement Learning for Dexterous Manipulation[END_REF], [START_REF] Ho | Generative adversarial imitation learning[END_REF], [START_REF] Kumar | Learning dexterous manipulation policies from experience and imitation[END_REF]. BC typically requires a large number of demonstrations and has issues with not observed trajectories. While these problems might be solved with additional expert supervision [START_REF] Ross | Reinforcement and imitation learning via interactive no-regret learning[END_REF] or noise injection in expert demonstrations [START_REF] Laskey | DART: Noise injection for robust imitation learning[END_REF], we address them by improving the standard BC framework. We use recent state-of-the-art CNN architectures and data augmentation for expert trajectories. This permits to significantly reduce the number of required demonstrations and to improve performance. Reinforcement learning (RL). RL learns to solve a tasks without demonstrations using exploration. Despite impressive results in several domains [START_REF] Mnih | Human-level control through deep reinforcement learning[END_REF], [START_REF] Silver | Mastering the game of Go with deep neural networks and tree search[END_REF], [17], [START_REF] Kober | Reinforcement learning in robotics: a survey[END_REF], RL methods show limited capabilities when operating in complex and sparse-reward environments common in robotics. Moreover, RL methods typically require prohibitively large amounts of interactions with the environment during training. Hierarchical RL (HRL) methods alleviate some of these problems by learning a high-level policy modulating low-level workers. HRL approaches are generally based either on options [START_REF] Sutton | Between mdps and semimdps: A framework for temporal abstraction in reinforcement learning[END_REF] or a feudal framework [27]. The option methods learn a master policy that switches between separate skill policies [START_REF] Frans | Meta learning shared hierarchies[END_REF]- [START_REF] Florensa | Stochastic neural networks for hierarchical reinforcement learning[END_REF]. The feudal approaches learn a master policy that modulates a low-level policy by a control signal [START_REF] Haarnoja | Latent space policies for hierarchical reinforcement learning[END_REF]- [START_REF] Hausman | Learning an embedding space for transferable robot skills[END_REF]. Our approach is based on options but in contrast to the cited methods, we pretrain the skills with IL. This allows us to solve complex and sparse reward problems using significantly less interactions with the environment during training. Combining RL and IL. A number of approaches combining RL and IL have been introduced recently. Gao et al. [START_REF] Gao | Reinforcement learning from imperfect demonstrations[END_REF] use demonstrations to initialize the RL agent. In [START_REF] Cheng | Fast policy learning through imitation and reinforcement[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF] RL is used to improve expert demonstrations, but does not learn hierarchical policies. Demonstrations have also been used to define RL objective functions [START_REF] Hester | Deep Q-Learning from demonstrations[END_REF], [START_REF] Nair | Overcoming exploration in reinforcement learning with demonstrations[END_REF] and rewards [START_REF] Zhu | Reinforcement and imitation learning for diverse visuomotor skills[END_REF]. Das et al. [START_REF] Das | Neural modular control for embodied question answering[END_REF] combine IL and RL to learn a hierarchical policy. Unlike our method, however, [START_REF] Das | Neural modular control for embodied question answering[END_REF] requires full task demonstrations and task-specific reward engineering. Moreover, the addressed navigation problem in [START_REF] Das | Neural modular control for embodied question answering[END_REF] has a much lower time horizon compared to our tasks. [START_REF] Das | Neural modular control for embodied question answering[END_REF] also relies on pre-trained CNN representations which limits its application domain. Le at al. [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF] train low-level skills with RL, while using demonstrations to switch between skills. In a reverse manner, we use IL to learn low-level control and then deploy RL to find appropriate sequences of pre-trained skills. The advantage is that our method can learn a variety of complex manipulations without full task demonstrations. Moreover, [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF] learn discrete actions and cannot be directly applied to robotics manipulations that require continous control.

In summary, none of the methods [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF], [START_REF] Cheng | Fast policy learning through imitation and reinforcement[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF] is directly suitable for learning complex robotic manipulations due to requirements of dense rewards [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF] and state inputs [START_REF] Cheng | Fast policy learning through imitation and reinforcement[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF], limitations to short horizons and discrete actions [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF], the requirement of full task demonstrations [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Le | Hierarchical imitation and reinforcement learning[END_REF], [START_REF] Cheng | Fast policy learning through imitation and reinforcement[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF] and the lack of learning of visual representations [START_REF] Das | Neural modular control for embodied question answering[END_REF], [START_REF] Cheng | Fast policy learning through imitation and reinforcement[END_REF], [START_REF] Sun | Truncated horizon policy search: Combining reinforcement learning & imitation learning[END_REF]. Moreover, our skills learned from synthetic demonstrated trajectories outperform RL based methods, see Section V-D.

III. APPROACH

Our RLBC approach aims to learn multi-step policies by combining reinforcement learning (RL) and pre-trained skills obtained with behavioral cloning (BC). We present BC and RLBC in Sections III-A and III-B. Implementation details are given in Section III-C.

A. Skill learning with behavioral cloning

Our first goal is to learn basic skills that can be composed into more complex policies. Given observation-action pairs D = {(o t , a t )} along expert trajectories, we follow the behavioral cloning approach [START_REF] Pomerleau | ALVINN: An autonomous land vehicle in a neural network[END_REF] and learn a function approximating the conditional distribution of the expert policy π E (a t |o t ) controlling a robot arm. Our observations o t ∈ O = R H×W×M are sequences of the last M depth frames. Actions a t = (v t , ω t , g t ), a t ∈ A BC are defined by the end-effector linear velocity v t ∈ R 3 and angular velocity ω t ∈ R 3 as well as the gripper openness state g t ∈ {0, 1}.

We learn the deterministic skill policies π s : O → A BC approximating the expert policy π E . Given observations o t with corresponding expert (ground truth) actions a t = (v t , ω t , g t ), we represent π s with a convolutional neural network (CNN) and learn network parameters (θ, η)

such that predicted actions π s (o t ) = (v t , ωt , ĝt ) minimize the loss L BC (π s (o t ), a t ) = λ [v t , ωt ] -[v t , ω t ] 2 2 + (1 - λ) g t log ĝt + (1 -g t ) log (1 -ĝt )
, where λ ∈ [0, 1] is a scaling factor which we empirically set to 0.9.

Our network architecture is presented in Figure 1(right). When training a skill policy π i s , such as reaching, grasping or pouring, we condition the network on the skill using the recent FiLM architecture [START_REF] Perez | Film: Visual reasoning with a general conditioning layer[END_REF]. Given the one-hot encoding s i of a skill i, we use s i as input to the FiLM generator which performs affine transformations of the network feature maps. FiLM conditions the network on performing a given skill, which permits learning a shared representation for all skills. Given an observation o t , the network CN N (θ|s i ) generates a feature map x i t conditioned on skill i. The spatially-averaged x i is linearly mapped with F C(η) to the action of π i s .

B. RLBC approach

We wish to solve composite manipulations without full expert demonstrations and with a single sparse reward. For this purpose we rely on a high-level master policy π m controlling the pre-trained skill policies π s at a coarse timescale. To learn π m , we follow the standard formulation of reinforcement learning and maximize the expected return E π ∞ k=0 γ k r t+k given rewards r t . Our reward function is sparse and returns 1 upon successful termination of the task and 0 otherwise. The RL master policy π m : O × A RL → [0, 1] chooses one of the K skill policies to execute the low-level control, i.e., the action space of π m is discrete: A RL = {1, . . . , K}. Note, that our sparse reward function makes the learning of deep visual representations challenging. We, therefore, train π m using visual features x i t obtained from the BC pretrained CN N (θ|s i ). Given an observation o t , we use the concatenation of skill-conditioned features {x 1 t , . . . , x K t } as input for the master CN N (µ), see Figure 1(right).

To solve composite tasks with sparse rewards, we use a coarse timescale for the master policy. The selected skill policy controls the robot for n consecutive time-steps before the master policy is activated again to choose a new skill. This allows the master to focus on high-level task planning rather than low-level motion planning achieved by the skills. We expect the master policy to recover from unexpected events, for example, if an object slips out of a gripper, by re-activating an appropriate skill policy. Our combination of the master and skill policies is illustrated in Figure 1(left).

RLBC algorithm. The pseudo-code for the proposed approach is shown in Algorithm 1. The algorithm can be divided into three main steps. First, we collect a dataset of expert trajectories D k for each skill policy π k s . For each policy, we use an expert script that has an access to the full state of the environment. Next, we train a set of skill policies {π 1 s , . . . , π K s }. We sample a batch of state-action pairs and update parameters of convolutional layers θ and the skills linear layer parameters η. Finally, we learn the master π m using the pretrained skill policies and the frozen parameters θ. We collect episode rollouts by first choosing a skill policy with the master and then applying the selected skill to the environment for n time-steps. We update the master policy weights µ to maximize the expected sum of rewards. C. Approach details Skill learning with BC. We use ResNet-18 for the CN N (θ|s i ), which we compare to VGG16 and ResNet-101 in Section V-A. We augment input depth frames with random translations, rotations and crops. We also perform viewpoint augmentation and sample the camera positions on a section of a sphere centered on the robot and with a radius of 1.40 m. We uniformly sample the yaw angle in [-15°, 15°], the pitch angle in [15°, 30°], and the distance to the robot base in [1.35m, 1.50m]. The impact of both augmentations is evaluated in Section V-B. We normalize the ground truth of the expert actions to have zero mean and a unit variance and normalize the depth values of input frames to [-1, 1]. We learn BC skills using Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with the learning rate 10 -3 and a batch size 64. We also use Batch Normalization [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF].

(o t ,a t )∈D k L BC (π k s (ot),
Task learning with RL. We learn the master policies with the PPO [START_REF] Schulman | Proximal policy optimization algorithms[END_REF] algorithm using the open-source implementation [START_REF] Kostrikov | Pytorch implementations of reinforcement learning algorithms[END_REF] where we set the entropy coefficient to 0.05, the value loss coefficient to 1, and use 8 episode rollouts for the PPO update. For the RLBC method, the concatenated skill features {x 1 t , . . . , x K t } are processed with the master network CN N (µ) having 2 convolutional layers with 64 filters of size 3 × 3. During pre-training of skill policies we update the parameters (θ, η). When training the master policy, we only update µ while keeping (θ, η) parameters fixed. We train RLBC using 8 different random seeds in parallel and evaluate the best one.

Real robot transfer. To deploy our method on the real robot, we use a state-of-the-art technique of learning sim2real transfer based on data augmentation with domain randomization [START_REF] Pashevich | Learning to augment synthetic images for sim2real policy transfer[END_REF]. This method uses a proxy task of cube position prediction and a set of basic image transformations to learn a sim2real data augmentation function for depth images. We augment the depth frames from synthetic expert demonstrations with this method and, then, train skill policies. Once the skill policy is trained on these augmented simulation images, it is directly used on the real robot.

IV. EXPERIMENTAL SETUP This section describes the setup used to evaluate our approach. First, we present the robot environment and the different tasks in Sections IV-A and IV-B. Next, we describe the synthetic dataset generation and skill definition for each task in Sections IV-C and IV-D.

A. Robot and agent environment

For our experiments we use a 6-DoF UR5 robotic arm with a 3 finger Robotiq gripper, see Figure 2. In simulation, we model the robot with the pybullet physics simulator [START_REF] Courmans | PyBullet, Python module for physics simulation, robotics and machine learning[END_REF]. For observation, we record depth images with the Microsoft Kinect 2 placed in front of the arm. The agent takes as input the three last depth frames o t ∈ R 224×224×3 and commands the robot with an action a t ∈ R 7 . The control is performed at 10 Hz frequency.

B. UR5 tasks

For evaluation, we consider 3 tasks: UR5-Pick, UR5-Bowl and UR5-Breakfast. The UR5-Pick task picks up a cube of a size between 3.5 cm and 8.0 cm and lifts it up, see Figure 2a. In UR5-Bowl the robot has to grasp the cube and place it in the bowl, see Figure 2b. The UR5-Breakfast task contains a cup, a bottle and a bowl as shown in Figure 2c. We use distinct ShapeNet [START_REF] Chang | ShapeNet: An information-rich 3D model repository[END_REF] object instances for the training and test sets (27 bottles, 17 cups, 32 bowls in each set). The robot needs to pour ingredients from the cup and the bottle in the bowl. In all tasks, the reward is positive if and only if the task goal is reached. The maximum episode lengths are 200, 600, and 2000 time-steps for UR5-Pick, UR5-Bowl, and UR5-Breakfast correspondingly.

C. Synthetic datasets

We use the simulated environments to create a synthetic training and test set. For all our experiments, we collect with random initial configurations where the objects and the end-effector are allocated within a workspace of 80 × 40 × 20 cm 3 . The synthetic demonstrations are collected using an expert script designed for each skill. The script has access to the full state of the system including the states of the robot and the objects. To generate synthetic demonstrations, we program end-effector trajectories and use inverse kinematics (IK) to generate corresponding trajectories in the robot joints space. Each demonstration consists of multiple pairs of the three last camera observations and the robot control command performed by the expert script. For UR5-Pick, we collect 1000 synthetic demonstrated trajectories for training. For UR5-Bowl and UR5-Breakfast, we collect a training dataset of 250 synthetic demonstrations. For evaluation of each task, we use 100 different initial configurations in simulation and 20 trials on the real robot. D. Skill definition UR5-Pick task is defined as a single skill. For UR5-Bowl and UR5-Breakfast, we consider a set of skills defined by expert scripts. For UR5-Bowl, we define four skills: (a) go to the cube, (b) go down and grasp, (c) go up, and (d) go to the bowl and open the gripper. For UR5-Breakfast, we define four skills: (a) go to the bottle, (b) go to the cup, (c) grasp an object and pour it to the bowl, and (d) release the held object. We emphasize that the expert dataset does not contain full task demonstrations and that all our training is done in simulation. When training the RL master, we execute selected skills for 60 consecutive time-steps for the UR5-Bowl task and 220 time-steps for the UR5-Breakfast task.

V. EVALUATION OF BC SKILL LEARNING

This section evaluates the different parameters of the BC skill training for the UR5-Pick task and a comparison with the state of the art. First, we evaluate the impact of the CNN architecture and data augmentation on the skill performance in Sections V-A and V-B. Then, we show that the learned policies transfer to a real-robot in Section V-C. Finally, we compare the BC skills with the state of the art in Section V-D.

A. CNN architecture for BC skill learning

Given the simulated UR5-Pick task illustrated in Figure 2a(left), we compare BC skill networks trained with different CNN architectures and varying number of expert demonstrations. Table I compares the success rates of policies with VGG and ResNet architectures. Policies based on the VGG architecture [START_REF] Simonyan | Very deep convolutional networks for large-scale image recognition[END_REF] obtain success rate below 40% with 100 training demonstrations and reach 95% with 1000 demonstrations. ResNet [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] based policies have a success rate above 60% when trained on a dataset of 100 demonstrations and reach 100% with 1000 demonstrations. Overall ResNet-101 has the best performance closely followed by ResNet-18 and outperforms VGG significantly. To conclude, we find that the network architecture has a fundamental impact on the BC performance. In the following experiments we use ResNet-18 as it presents a good trade-off between performance and training time.

When examining why VGG-based BC has a lower success rate, we observe that it has higher validation errors compared to ResNet. This indicates that VGG performs worse on the level of individual steps and is hence expected to result in higher compounding errors.

B. Evaluation of data augmentation

We evaluate the impact of different types of data augmentations in Table II. We compare training without data augmentation with 3 variants: (1) random translations, rotations and crops, as is standard for object detection, (2) record each expert synthetic demonstration from 10 varying viewpoints and (3) the combination of ( 1) and (2).

Success rates for UR5-Pick on datasets with 20, 50 and 100 demonstrations are reported in Table II. We observe that data augmentation is particularly important when only a few demonstrations are available. For 20 demonstrations, the policy trained with no augmentation performs at 1% while the policy trained with standard and viewpoint augmentations together performs at 75%. The policy trained with a combination of both augmentation types performs the best and achieves 93% and 100% success rate for 50 and 100 demonstrations respectively. In summary, data augmentation allows a significant reduction in the number of expert trajectories required to solve the task.

C. Real robot experiments

We evaluate our method on the real-world UR5-Pick illustrated in Figure 2a(right). We collect demonstrated trajectories in simulation and train the BC skills network applying standard, viewpoint and sim2real augmentations. We show that our approach transfers well to the real robot using no real images. The learned policy manages to pick up cubes of 3 different sizes correctly in 20 out of 20 trials.

D. Comparison with state-of-the-art methods

One of the few test-beds for robotic manipulation is FetchPickPlace from OpenAI Gym [START_REF] Plappert | Multi-goal reinforcement learning: Challenging robotics environments and request for research[END_REF] implemented in mujoco [START_REF] Todorov | MuJoCo: A physics engine for model-based control[END_REF], see Figure 3a. The goal for the agent is to pick up the cube and to move it to the red target (see Figure 3a). The agent observes the three last RGB-D images from a camera placed in front of the robot o t ∈ R 100×100×4×3 . The positions of the cube and the target are set at random for each trial. The reward of the task is a single sparse reward of success. The maximum length of the task is 50 time-steps.

For a fair comparison with [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF], we do not use any data augmentation. We report the success rate of ResNet-18 policy in Figure 3b. We follow [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF] and plot the success rate of both RL and IL methods with respect to the number of episodes used (either trial episodes or demonstrations). Our approach outperforms the policies trained with an imitation learning method DAgger [START_REF] Ross | Reinforcement and imitation learning via interactive no-regret learning[END_REF] in terms of performance and RL methods such as HER [START_REF] Andrychowicz | Hindsight experience replay[END_REF] and DDPG [START_REF] Lillicrap | Continuous control with deep reinforcement learning[END_REF] in terms of data-efficiency. According to [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF], DAgger does not reach 100% even after 8 * 10 4 demonstrations despite the fact that it requires an expert during training. HER reaches the success rate of 100% but requires about 4 * 10 4 trial episodes. Our approach achieves the 96% success rate using 10 4 demonstrations.

Our policies differ from [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF] mainly in the CNN architecture. Pinto et al. [START_REF] Pinto | Asymmetric Actor Critic for Image-Based Robot Learning[END_REF] use a simple CNN with 4 convolutional layers while we use ResNet-18. Results of this section confirm the large impact of the CNN architecture on the performance of visual BC policies, as was already observed in Table I.

VI. EVALUATION OF RLBC

This section evaluates the proposed RLBC approach and compares it to baselines introduced in Section VI-A. First, we evaluate our method on UR5-Bowl in Section IV-B. We then test the robustness of our approach to various perturbations such as dynamic changes of object positions, dynamic occlusions, unseen object instances and the increased probability of collisions due to small distances between objects. We show that RLBC outperforms the baselines on those scenarios both in simulation and on a real robot. Note, that our real robot experiments are performed with skills and master policies that have been trained exclusively in simulation using sim2real augmentation [START_REF] Pashevich | Learning to augment synthetic images for sim2real policy transfer[END_REF]. We use the same policies for all perturbation scenarios. Qualitative results of our method are available in the supplementary video.

A. Baseline methods

We compare RLBC with 3 baselines: (a) a fixed sequence of BC skills following the manually pre-defined correct order (BC-ordered); (b) an open-loop controller estimating positions of objects and executing an expert script (Detect & Plan); (c) a closed-loop controller performing the same estimationbased control and replanning in case if object positions change substantially (Detect & Replan). We use the same set of skills for RLBC and BC-ordered. We train the position estimation network using a dataset of 20.000 synthetic depth images with randomized object positions. All networks use ResNet-18 architecture and are trained with the standard, viewpoint and sim2real augmentations described in Section V-B.

UR5-Bowl perturbations

B. Results on UR5-Bowl with no perturbations

We first evaluate RLBC and the three baselines on the UR5-Bowl task (see Figure 2b). When tested in simulation, all the baselines and RLBC manage to perfectly solve the task. On the real-world UR5-Bowl task, BC-ordered and Detect & Plan baselines sometimes fail to grasp the object which leads to task failures (see Table III, first row). On the contrary, RLBC solves the task in all 20 episodes given its ability to re-plan the task in the cases of failed skills.

We have also attempted to solve the simulated UR5-Bowl task without skills by learning an RL policy performing lowlevel control. We have used ImageNet pre-trained ResNet-18 to generate visual features. The features were then used to train low-level RL control policy with PPO. Whereas such a low-level RL policy did not solve the task a single time after 10 4 episodes, RLBC reaches 100% after 400 episodes.

C. Robustness to perturbations

Robustness to dynamic changes in object position. We evaluate RLBC against the baselines in the UR5-Bowl scenario where the cube is moved several times during the episode. We plot success rates evaluated in simulation with respect to the number of position changes in Figure 4a. We observe the stability of RLBC and the fast degradation of all baselines. As both RLBC and BC-ordered use the same set of skills, the stability of RLBC comes from the learned skill combination. The "Moving objects" row in Table III reports results for 3 moves of the cube evaluated on the real robot. Similarly to the simulated results, we observe excellent results of RLBC and the degraded performance for all the baselines. Robustness to occlusions. We evaluate the success of UR5-Bowl task under occlusions. Each occlusion lasts 3 seconds and covers a large random part of the workspace by a cardboard. Figure 4b shows success rates with respect to the number of occlusions in the simulated UR5-Bowl environment. Similarly to the perturbation results in Figure 4a, RLBC demonstrates high robustness to occlusions while the performance of other methods quickly degrades. The "Occlusions" row in Table III reports results for a single occlusion performed during the real-robot evaluation. Baseline methods are strongly influenced by occlusions except Detect & Plan which performs well unless occlusion happens during the first frames. Our RLBC policy performs best compared to other methods. Robustness to new object instances. We evaluate the robustness of methods to the substitution of a cube by other objects not seen during the training of UR5-Bowl task. Figure 4c shows the success rate of RLBC and other methods with respect to the number of new objects in simulation. The novel objects are ordered by their dissimilarity with the cube. The difficulty of grasping unseen objects degrades the performance of grasping skills. In contrast to other methods RLBC is able to automatically recover from errors by making several grasping attempts. Table III reports corresponding results on a real robot where the cube has been replaced by 10 unseen objects. Similarly to the other perturbations we observe superior performance of RLBC. Impact of the distance between objects. We vary the distance between a bottle and a cup in the UR5-Breakfast task. The smaller distance between objects A and B implies higher probability of collision between a robot and A when grasping B behind A. The choice of the grasping order becomes important in such situations. While our method is able to learn the appropriate grasping order to maximize the chance of completing the task, the BC-ordered and other baselines use pre-defined order. Figure 4d demonstrates the performance of RLBC and BC-ordered for different object distances in the simulated UR5-Breakfast task. As expected, RLBC learns the correct grasping order and avoids most of collisions. The performance of BC-ordered strongly degrades with the decreasing distance. In the real-world evaluation, both RLBC and ordered skills succeed in 16 out of 20 episodes when the distance between objects is larger than 10 cm. However, the performance of BC-ordered drops to 8/20 when the cup and the bottle are at 4cm from each other. In contrast, RLBC chooses the appropriate object to avoid collisions and succeeds in 16 out of 20 trials.

VII. CONCLUSION

This paper presents a method to learn combinations of primitive skills. Our method requires no full-task demonstrations nor intermediate rewards and shows excellent results in simulation and on a real robot. We demonstrate the versatility of our approach in challenging settings with dynamic scene changes. Future work will include learning multiple tasks with shared skills and addressing contact-rich tasks.

Fig. 1 :

 1 Fig. 1: Illustration of our approach. (Left): Temporal hierarchy of master and skill policies. The master policy πm is executed at a coarse interval of n time-steps to select among K skill policies π 1 s . . . π K s . Each skill policy generates control for a primitive action such as grasping or pouring. (Right): CNN architecture used for the skill and master policies.

  ** Collect expert data *** 2: for k ∈ {1, . . . , K} do 3: Collect an expert dataset D k for the skill policy π k s 4: *** Train {π 1 s , . . . , π K s } by solving: *** 5: θ, η = arg min θ,η K k=1

2 :

 2 UR5 tasks used for evaluation: (a) task of picking up the cube, (b) task of bringing the cube to the bowl, (c) task of pouring the cup and the bottle into the bowl. (Left) simulation, (right) real robot. (a) FetchPickPlace (b) Comparison with [7] Fig. 3: Comparison of BC ResNet-18 with state of the art [7] on the FetchPickPlace task. BC ResNet-18 results are reported for 200 different initial configurations.

Fig. 4 :

 4 Fig. 4: Performance of RLBC and baseline methods in simulated environments under perturbations: (a) Dynamic changes of cube position; (b) Dynamic occlusions; (c) Replacing the cube by unseen objects; (d) Decreasing the distance between objects.

TABLE I :

 I Evaluation of BC skills trained with different CNN architectures and number of demonstrations on the UR5-Pick task in simulation.

	Demos None Standard Viewpoint	Standard & Viewpoint
	20	1%	49%	39%	75%
	50	5%	81%	79%	93%
	100	65%	97%	100%	100%

TABLE II :

 II Evaluation of ResNet-18 BC skills trained with different data augmentations on UR5-Pick task in simulation.

TABLE III :

 III Comparison of RLBC with 3 baselines on the real-world UR5-Bowl task with dynamic changes of the cube position, dynamic occlusions and new object instances.

VIII. ACKNOWLEDGMENTS This work was partially supported by the HPC resources from GENCI-IDRIS (Grant 20XX-AD011011163), Louis Vuitton ENS Chair on Artificial Intelligence, and the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).


E-M Abbou-Ou-Cherif

E Piron

A Chateauneuf

D Miclet

S Villette
email: sylvain.villette@agrosupdijon.fr

On-the-field simulation of fertilizer spreading: Part 3 -Control of disk inclination for uniform application on undulating fields

To improve the quality of fertiliser centrifugal spreading, several control devices have already been developed to manage some disruptions occurring on horizontal fields. However, controls which manage disruptions on non-flat fields due to changes in tractor attitude (pitch and roll) have not been developed. In this study, the design of a new control device is developed for a twin-disk spreader by considering two new degrees of freedom for each disk and controlling the longitudinal and lateral tilts of each disk. These tilt corrections are derived from solving a constrained optimization problem. The cost function is the weighted sum of squared differences between the travelled distances of the particles obtained in the considered topography and those that would have been obtained on a horizontal surface. The weighting coefficients are provided by the experimental angular mass flow distribution. In order to reduce the computation time and expect the use of the method for a real-time correction device, a simplified optimization problem is proposed. The method is assessed for three ground surface con-figurations: longitudinal slope break, side slope break and combined slopes. The simulation demonstrates that using the new control device the application errors are lower than 10% while they reach up to 40% without tilt correction. Moreover the study shows that the problem of optimization uniformity on non-flat fields can be solved computing the travelled distance for only a very few number of particles of mean diameter. This helps in reducing the computational time required to solve the optimization problem and in making possible the de-velopment of real-time control devices.

Introduction

In agriculture, centrifugal spreaders are widely used for fertilizer application, with up to 90% of the market share in Europe [START_REF] Van Liedekerke | DEM simulations of the particle flow on a centrifugal fertilizer spreader[END_REF]. In the last decades, numerous sensors and control mechanisms have been developed and marketed to improve spreading accuracy and environmental compliance [START_REF] Yule | New spreading technologies for improved accuracy and environmental compliance[END_REF]. Various new technologies have been used such as on-board weighing, variable rate application device, border management device, variable width spreading device, headland management system, start and stop control device, or self-calibration system. For instance, modern centrifugal spreaders can adjust the fertilization rate using variable rate technologies (VRT), based on field information provided by geographical information systems (GIS). Drop point and flow rate control devices are also commonly used to adapt the shape of the spread pattern deposition to the two-dimensional geometry of the field. Using simulations, [START_REF] Virin | Application of optimization techniques for an optimal fertilization by centrifugal spreading[END_REF] already demonstrated that such control systems improve the accuracy of the spreading, keeping the applied rate in a range of +/-10% of the local targeted application for flat fields. Nevertheless, when spreaders are used on non-horizontal surfaces, the accuracy of the spreading may be called into question.

Few works addressed the effect of spreader level on the fertilizer distribution uniformity. [START_REF] Parish | Effect of impeller angle on pattern uniformity[END_REF] studied the effect of the longitudinal tilt (front to rear angle) of the spreader and demonstrated that small changes such as 5 degrees caused significant pattern distortion. [START_REF] Yildirim | Effect of side to side spreader angle on pattern uniformity in singleand twin-disc rotary fertilizer spreaders[END_REF] studied the effect of the transverse tilt (side to side angle) and demonstrated that the fertilizer distribution pattern was skewed when the spreader was not horizontal. These studies highlighted that an automatic leveling system would be required to maintain an even fertilizer distribution when the ground surface was horizontal but they did not consider the spreading distortions in the case of non-flat fields. More recently, [START_REF] Grafton | The effect of field conditions on in-field spread patterns from twin disk spreaders[END_REF] pointed out the negative effect of slope on the spreading quality because of a skewing of the spread pattern. [START_REF] Horrocks | Developing guidelines for fertiliser spreading on west coast humps and hollows[END_REF] also demonstrated that spreading fertilizer on humps and hollows affects the distribution patterns and the uniformity of fertilizer application depending on the bout width, and the driving position with respect to humps, slopes and hollows.

Considering the proportion of the arable land located on slopes and hilly regions, the proportion of field areas that could be affected by nonuniform fertilizer application cannot be ignored. For example, in the hilly region of Bend Subcarpathians, in Romania, [START_REF] Chendes | Analyzing the landforms-agricultural land-use types relationship using a DTM-based indicator[END_REF] analysed that 52% of arable lands are situated on low slope lands and 18% flat slow lands.

Thus, in order to improve the quality of the spreading, a new active control device needs to be developed to ensure even fertilizer distribution on undulating fields regardless of their slope. To develop this new technical solution, a traditional approach would lead to costly or unrealistic experiments because of the diversity of slope situations and the difficulty of in-field measurements. Thus, virtual experiments based on numerical simulations were preferred to overcome the limits of infield experiments. The use of hybrid approaches [START_REF] Reumers | Experimental characterisation of the cylindrical distribution pattern of centrifugal fertiliser spreaders: towards an alternative for spreading measurements[END_REF][START_REF] Villette | Hybrid centrifugal spreading model to study the fertiliser spatial distribution and its assessment using the transverse coefficient of variation[END_REF] was of particular interest to combine some theoretical motion models with experimental data measured at some steps of the spreading process. In order to obtain realistic results, [START_REF] Villette | Hybrid centrifugal spreading model to study the fertiliser spatial distribution and its assessment using the transverse coefficient of variation[END_REF] suggested the use of the Hybrid Centrifugal Spreading Model which combined a mechanistic approach based on the use of mechanical relationships and a stochastic approach based on the use of statistical distribution input parameters (derived from actual experiments). Nevertheless, the model described in [START_REF] Villette | Hybrid centrifugal spreading model to study the fertiliser spatial distribution and its assessment using the transverse coefficient of variation[END_REF] was developed for flat ground conditions and did not take into account any disruptions due to the slope in the case of uneven or sloped ground (i.e. non-flat fields).

In order to study the effect of the slope on the quality of fertilizer spreading and design a correcting active control device, three main steps were required.

First, a new model was developed (Abbou-ou-cherif et al., 2017b) to simulate fertilizer spreading during the motion of a tractor on non-flat fields simulated by a digital elevation model. The coordinate system, the geometrical parametrization of tractor and spreader, and the ballistic flight model required for this study had been described by Abbouou-cherif et al. (2017b). In this first step, various experimental measurements were also carried out to provide the input data required by the model for ammonium nitrate. The particle size distribution was derived from a sieve analysis conducted in accordance with the standard ISO 8397 (1988). The horizontal outlet angle and vertical mass flow distribution was deduced from previous measurements [START_REF] Villette | Estimation of two-dimensional fertiliser mass flow distributions by recording granule impacts[END_REF]. The horizontal angular mass flow distribution and the particle drag coefficient were deduced from 2D spread patterns measured with a rotating test bench [START_REF] Piron | CEMIB: an innovative bench for spreader ecodesign[END_REF] for a concave disk equipped with two radial vanes (radius: 395 mm). The model was then assessed and validated by comparing simulated and measured 2D spread patterns for three disk rotation speeds (600, 800 and 1000 rpm) and three disk tilt angles (0°, 5°, 10°).

Second, the model was used (Abbou-ou-cherif et al., 2017a) to compute the local application rate for various simulated 3D slope configurations, to quantify the resulting over-and under-application, and to understand the mechanisms of spread pattern distortions. The use of the hybrid model ensured the study of application uniformity in theoretical but realistic conditions. Various situations had been investigated such as regular longitudinal or side slope and irregular slope with longitudinal or side slope break (i.e. change of slope). Application rate maps were simulated and over-and under-application areas were characterized by computing the mean longitudinal or transverse application rates as follows:

= × Q 100 mean m N M N i j N ij T 1 (1)
where i is the index of a virtual straight row of collecting trays of 1 × 1 m (in the longitudinal or transverse direction), m ij is the accu- mulated mass of particles in the jth collecting tray of the ith row, N is the constant number of collecting trays in every row, M is the total mass of particles in the considered area, N T is the total number of collecting trays in the same area and Q meani is the mean application rate in the ith row.

In the case of regular non-flat fields (with constant slope), regarding the mean transverse application rates, Abbou-ou-cherif et al. (2017a) computed application errors within ± 10% even at steep slopes. However, in the case of a longitudinal change of slope (for 10°positive slope) the deviation from the expected uniform rate reached a maximum of +45% and a minimum of -25% (deviation observed on the mean longitudinal application rate). These deviations from the targeted rate clearly demonstrated the need for correcting systems to keep the application variation within an acceptable range. Considering the level of the spreader, when working on a horizontal surface, [START_REF] Parish | Effect of impeller angle on pattern uniformity[END_REF] and [START_REF] Yildirim | Effect of side to side spreader angle on pattern uniformity in singleand twin-disc rotary fertilizer spreaders[END_REF] already suggested the design of an automatic levelling system to ensure an even fertilizer distribution. In the case of non-horizontal fields, the simulation results obtained by Abbou-oucherif et al. (2017a) not only help in understanding why the spread pattern is affected by the slope (in terms of tractor attitude and in terms of intersection between the ground surface and particle trajectories) but also provide help in designing appropriate correcting devices to prevent application errors.

Third, the model is now used to design a new correction system. Thus, the objective of this paper is to design an active control device for dual-disk centrifugal spreader to dynamically adapt the disk tilt to the site-specific slope and thus limit application errors when fertilizers are spread onto non-flat fields, regardless of the tractor path and ground slope configuration. The main contributions of the paper are: (1) the design of an active system to control the longitudinal and the lateral tilt of each spinning disk independently, (2) the design of a control strategy derived from an optimization problem which is adapted to reduce the computational time and consider the development of real-time control, (3) the assessment of the control strategy using numerical simulations. -ou-cherif et al. (2017a) demonstrated that fertilizer application errors were particularly significant in the case of irregular nonflat fields (i.e. with longitudinal or transverse change of slope). This is due to the difference between the inclination of the tractor and the ground surface that leads to important deviations in the lengths of the particle ballistic flights. These deviations also take place on regular fields with longitudinal or side slope but to a less extent. For a given situation (e.g. a particular inclination of the tractor and of the ground surface) the value of these deviations is not the same in all the directions and depends on the angular direction of the particle around the disk (see Abbou-ou-cherif et al. (2017a) for details). These changes in the lengths of the particle ballistic flights directly affect the shape of the spread pattern deposition and, consequently, affect the application uniformity. Thus, it appears that the control of the length of the ballistic flights is required to achieve a return to a spread pattern deposition close to the one that would be obtained on a horizontal surface and thus to keep the application variation within an acceptable range. Therefore, a method of correcting application errors on non-flat field surfaces is developed in this paper by controlling the lengths of the ballistic flights of the particles during the spreading process.

Formulation of the optimization problem

Abbou

Optimal control strategy

Changing particle flight lengths with current marketed centrifugal spreaders can be achieved through several means such as varying the disk velocity, the disk cone angle or the length of vanes. However, these parameters affect the flight lengths of all the particles in the whole spread pattern, and is not adapted to expect a differentiated correction depending on the angular direction of the particles around the spinning disk. As a consequence, a new control system has been developed in this study. It is based on the introduction of two new degrees of freedom for each disk (Fig. 1) so that its orientation can be varied in a three-dimensional space and specified in a spherical coordinate system. Thus, this orientation system can not only control the flight length of the particles, but also target the areas of application errors (i.e. a particular angular sector of spreading). Based on the system represented in Fig. 1, two control variables have been derived for each disk at each new position of the tractor on its trajectory: T (longitudinal tilt) and S (side or lateral tilt). For this purpose, a control strategy was implemented using the simulation model detailed in Abbou-ou-cherif et al. (2017a) and Abbou-ou-cherif et al. (2017b). This model computes the fertilizer rate applied on the field based on input data pertaining to the centrifugal spreader parameters, the particle characteristics and the digital elevation model (DEM) of the field. A classical control strategy would consist of inverting this direct model in order to derive the values of the control variables from the required application rate. However, this task is difficult to achieve because the intersections of the particle ballistic flights with the DEM can only be computed numerically. As a consequence, the control strategy chosen in this study consists in deriving the control variables by solving an optimization problem where the objective of the cost function is to minimize the sum of squared deviations of particle travelled distances compared to distances that would be obtained on a flat field. The control strategy developed in the following for one disk is applied separately for both disks.

The purpose of a combined action of longitudinal and side (i.e lateral) tilt is to correct the changes of the particle travelled distances due to ground topography, with respect to the distances that would be obtained on a horizontal surface. As illustrated in Fig. 2, the travelled distance is defined by the curvilinear length that separates the disk center from the landing position of a particle on the DEM (Abbou-oucherif et al., 2017a). Considering a particle with a travelled distance p 0i in a flat configuration, the same particle has a range p i on a non-flat field such as there is a non-zero error defined by:

= p p i i i 0 (2)
Then, the goal of the optimisation problem is to minimize the squared error i 2 by controlling the disk longitudinal and side tilts on non-flat fields.

Considering the results presented in Abbou-ou-cherif et al. (2017a) and their analysis, it has been assumed that the application errors could be reduced by keeping particle travelled distances equal to those obtained on a flat field, regardless of the tractor orientation. This objective can be mathematically expressed by a cost function which needs to be minimized at each point of the tractor trajectory. The cost function is as follows:

= = F T S d v T S ( , ) ( , , , , , ) c k N k k k k 1 2 (3)
where N is the number of particles ejected by the spinning disk at the considered tractor location, k is the angular direction of the kth particle measured with respect to the opposite of the travel direction, d k , v k and k are respectively the diameter, the outlet velocity and the vertical outlet angle (with respect to the disk plane) and the of the kth particle. The optimization output is a vector = x T S ( , )with the optimal va- lues of T and S to be used for the disk inclination. In practice, the value of i is in the interval [-120°, 120°] where a single spread pattern can be located. These limits have been deduced from numerous experimental measurements of single spread patterns captured with the rotating test bench CEMIB [START_REF] Piron | Centrifugal fertiliser spreaders: a new method for their evaluation and testing[END_REF] used in Irstea (the National Research Institute of Science and Technology for Environment and Agriculture, France) to design and certify centrifugal spreaders.

Spreading configuration

Typical values have been chosen for the spreading parameters. The simulations were computed for a tractor speed of 12 km/h, a working width of 28 m and a flow rate of 130 kg/ha. The spinning disks were equipped of radial vanes with a radius of 395 mm and a vertical angle of 13.25°. The disk rotation speed was 800 rpm. Considering particle characteristics, the density was 1750 kg m -3 , the drag coefficient was 0.54 and the mean diameter was 3.1 mm. The vertical and horizontal angular distributions of the flow rate around the disk were derived from experimental measurements (Abbou-ou-cherif et al., 2017b). These distributions were measured when the disk axle was vertical. Thus, in this work, it was assumed that these distributions (measured with respect to the disk) are still the same irrespective of the tilt of the disk.

Spatial discretization step

The numerical resolution of the optimization problem requires discretization in space along the tractor trajectory. With centrifugal spreaders, the fertilizer particles are discontinuously thrown from the disk. At each vane revolution, a set of particles is taken up by the vane, then undergoes a centrifugal acceleration and finally is thrown into the field when the particles reach the extremity of the vane. Thus, particles are periodically ejected from the spinning disk at each vane revolution. In order to take into account this discontinuous process, the coordinates of the particle landing points are computed after each disk revolution, at each position of the tractor along its trajectory. Considering the distance L crossed by the tractor, the distance L travelled by the tractor during the time required for one revolution of the spinning disk is as follows:

= L v 2 t d (4)
where v t (expressed in m s -1 ) is the tractor velocity and d is the ro- tational speed of the disk (expressed in rad s -1 ).

Considering a tractor velocity of 12 km/h and a disk rotation speed of 800 rpm, the distance L travelled by the tractor during one disk revolution is 0.25 m. This value is used in this study to discretize the tractor trajectory and thus define the trajectory points where the spread pattern needs to be computed.

Simplification of the cost function for real time computation

Deriving the optimal vector x * for each disk has to be done in real time at each point of the trajectory. However, the distance step L = 0.25 m is travelled in less than 0.075 s which is a time span significantly lower than the time necessary to simulate a single spread pattern, with at least 1000 particles (considering the usual working width, flow rate and velocity speed). As a consequence, the cost function in Eq. (3) needs to be simplified. Furthermore, to take into account the mechanical and dynamic limits of the actuators, the optimization problem has to be constrained.

The distortion of a single spread pattern can be represented by the displacements between the landing positions of the particles when they are spread on a horizontal surface and their landing positions on the non-flat field. For example, these displacements are illustrated in Fig. 3 in the case of a side slope break. Considering this example, two observations can be made. First, the travelled distance variation depends on the angular sector of spreading. In the case illustrated in Fig. 3a, the magnitude of the displacement vectors is zero for the particles ejected on the left side or at the rear of the tractor but the displacement can reach 4 m for the particles ejected on the right side. Moreover, the value of the displacement magnitude changes very progressively as a function of the angular direction (e.g. in Fig. 3, the displacement magnitude increases from the rear to the right side). Second, similar displacement fields are obtained whether by considering all the particles of the single spread pattern (Fig. 3a), or considering only a selection of particles (Fig. 3b) chosen with only one diameter (equals to the mean diameter: 3.1 mm) and thrown out with a constant outlet velocity (i.e. equals to the mean value: 44.2 m s -1 ) and a constant vertical outlet angle with respect to the disk plane (i.e. equals to the mean value: 8.7°).

Considering that the whole angular spread sector can be decomposed into a set of adjacent angular sectors of same angular width , the cost function described by Eq. ( 3) can be rewritten as follows: where n is the total number of angular sectors, N j is the number of particles to be considered in the jth sector.

= + + + = = = F T S d v T S d v T S d v T S ( , ) ( , , , , , ) ( , , , , ) . .. ( , , , , ) c k
N k k k k k N k k k k k N nk nk nk nk 1 2 1 1 1 1 1 2 2 2 2 2 1 2 n 1 2 (5)
Considering that the angular directions of all the particles are included in the interval [-120, 120°] (see Section 2.1), the value of jk is bounded as follows:

+ × < + × j j 120 ( 1) 120 jk (7)
which can also be written:

< + j jk j 1 (8)
with:

= + × j 120 ( 1) j (9)
with j from 1 to n.

As previously analysed (Fig. 3), the displacement computed for the whole range of particle characteristics can be estimated to the one computed with the mean characteristics of the particle (i.e. mean diameter and mean ejection parameters), the cost function can be simplified as follows:

= = = F T S d v T S ( , ) ( , ¯, ¯, ¯, , ) c j n k Nj jk 1 1 2 (10)
where d ¯, v ¯and ¯are the mean values of the particle diameter, the outlet velocity and the vertical outlet angle (with respect to the disk plane) respectively. Moreover, the previous analysis (Fig. 3) shows that the displacements can be considered as constant within one angular sector. Thus Eq. ( 10) yields:

= × = F T S Nj d v T S ( , ) ( ( , ¯, ¯, ¯, , )) c j n j 1 2 (11)
In order to express the cost function with respect to the proportion of fertiliser amount in each angular sector, a new expression of the cost function is derived as follows:

= = × = F T S F T S N Nj N d v T S ( , ) ( , ) ( , ¯, ¯, ¯, , ) c c j n jk 1 2 N (12) = = x x F ( ) ( ) c j n j j 1 2 N (13)
where x is the vector (T, S), x ( ) j 2 is computed for one particle considering the mean values of the particle diameter, the outlet velocity and the vertical outlet angle:

= x d v T S ( ) ( , ¯, ¯, ¯, , ) j jk 2 2 (14)
and the weighted coefficient j expresses the proportion of the fertilizer rate in the jth angular sector:

= N N j j (15)
Thus, in Eq. ( 13), the squared error obtained for the jth angular sector is weighted by j according to the fertilizer rate in the jth angular sector. Therefore, the weight of each angular sector in the single spread pattern distortion is taken into account. Using the density function of the horizontal distribution of particles around the disk, the coefficient j can be calculated as the probability for a particle an angular direction φ to be included in a given angular sector. Using the notation defined in Eq. ( 9), j is defined as follows:

= < = = + + P F F F f d ({ }) ( ) ( ) ( ) ( ) j i j a j a j a a 1 1 120 (16)
where

< + P ({ }) i j 1
is the probability of having in the interval ]φ i ,φ i+1 ], F is the cumulative distribution function and f a the probability density function of the particle horizontal angle.

In practice, j is directly derived from the experimental measurement of the horizontal angular distribution which is used as input data for the model.

Regarding the simplified cost function defined in Eq. ( 13), the higher the number n of angular sectors is, the closer the optimal vector x * becomes to the value that would be obtained using the cost function defined in Eq. ( 3). However, as illustrated in Fig. 3b, the trend of the displacements is a gradual and continuous change with respect to the angular direction. As a consequence, the choice of a small number of angular sectors and thus a small number of particles can be sufficient to carry out the optimization at each point of the tractor trajectory. In practice, the choice of 3 angular sectors is sufficient to provide good results as shown in Section 3.

Optimization constraints

To take into account the mechanical and dynamic limits of the actuators, several constraints have to be expressed as lower and upper bounds concerning the variables that have to be optimized (i.e. T and S) and their time derivatives. The set A of constraints is defined as follows:

L T U C L S U C L U C L U C ( 1) ( 2)
( 3) (17

)
where L bi and U bi are respectively the lower and upper boundary defining the constraint Ci.

As the trajectory of the tractor is discretized at every distance step L (Eq. ( 4)), the constraints C3 and C4 are rewritten to be discretized by taking into account the tractor velocity v t . For example, the constraint C3 is modified as follows:

L v dT dt v U v . 1 b t t b t 3 3 (18) L v dT dt dt dl U v . b t b t 3 3 (19)
where l and t are respectively the tractor displacement and the time in the continuous domain. Thus, in the discrete domain, the following relationships are established:

L v T T L U v b t k k b t 3 1 3 (20) + + L v L T T U v L T b t k k b t k 3 1 3 1 ( 21
)
where k is the index of the tractor position, T k-1 and T k are two successive values of T regarding the spatial discretisation.

Similarly, considering the constraint C4, bounds are established for the variable S at the location k as a function of its value found in the previous position k-1. The final formulation of the constrained optimization problem is given by:

= = + + + + = x x A L T U C L S U C L T T L T C L S S L S C min ( ) ( 1)
( 2)

( 3) ( 4)

x k A j n j j b k b b k b L v k k U v k L v k k U v k 1 2 1 1 2 2 1 1 1 1 b t b t b t b t 3 3 4 4 (22)
For numerical computation, the tilt angles of the actuators were limited at ± 30°(L b1 = L b2 = -30°, U b1 = U b2 = +30°) and the rotational speeds were limited at ± 4°/s (L b3 = L b4 = -4°/s, U b3 = U b4 = +4°/s).

Optimization problem solving and algorithm

The optimization problem was solved using the function fmincon of the optimization toolbox of the software Matlab (2005). This function allows finding a local optimum thanks to several evaluations of the cost function using the simulation model of the ballistic flights. In order to improve the computation time, an analytical expression of the cost function is derived at each new position of the trajectory. For this purpose, it is proposed to fit a third-order polynomial surface to a fixed number of points of the 2D cost function at the trajectory index k. These points are in the vicinity of the optimal solution x k 1 at the previous position of the trajectory. The polynomial surface coefficients are obtained by the method of least squares to solve the following linear system:

= S T S T S T S T T S S T a

b c d e f g h i j F T S [1 ] [ ( , )] p p p p p p p p p p p p c p p 2 2 2 2 3 3 N ( 23
)
where p is the pth point of fit (with p from 1 to n p the number of points used to locally fit the polynomial surface). In this work, n p was set at 10 points.

In order to calculate the cost function necessary to derive x * (k + 1), the points required to fit the cost function were then translated to the of the previously found solution x * Once the analytical expression of the polynomial surface was known, its gradient and Hessian matrix were computed and provided as additional arguments to the function fmincon. This approach ensured finding systematically a solution for the optimization problem at each point of the tractor trajectory and in a limited time frame. Moreover, the convexity of the cost function ensured that each solution was the optimal solution of the optimization problem.

Results and discussions

To assess the effectiveness of the control strategy based on the disk inclination, simulations were carried out for three ground surface configurations (i.e. longitudinal slope break, side slope break and undulating surface with combined slopes) using three discrete elevation models. These are presented in Fig. 4 (upper row). The investigated configurations are representative of the elevation changes of agricultural fields regarding the longitudinal and transverse slope with respect to the tractor travel direction. In each case, the fertilizer application maps were computed with (Fig. 4, lower row) and without (Fig. 4, middle row) the disk inclination control strategy. Considering longitudinal and side slope breaks, the digital elevation models used a horizontal surface and a 10°tilted surface, connected by a surface with a radius of curvature of 100 m. The undulating elevation model (with combined slopes) simulated undulating hills and valleys with slopes set between -6°and +6°. All these parameter values ensured the realistic nature of the simulations as explained in Abbou-ou-cherif et al. (2017a). When the control strategy is used, the application map presented in Fig. 4 (lower row) results from the use of only 3 angular sectors (n = 3) in the optimization problem. For the three discrete elevation models, comparing the middle and the lower row of Fig. 4 shows that, over-application and under-applications are much lower when the disk inclination control system is used. Thus, the corrections of the disk tilts significantly improve the uniformity of the fertilizer application for the three topographic reliefs simulated.

Fig. 5 presents the longitudinal and transverse application rates computed using Eq. ( 1) considering the areas delineated in Fig. 4 for the three ground surface configurations. When the control strategy was used, the mean application rates were computed considering various number n of angular sectors for the optimization problem. In all cases, Fig. 5 shows that the application error magnitudes are significantly reduced since they are within the interval ± 10% when the control strategy is used. Furthermore, Fig. 5 shows that choosing n = 3 provides good results and no uniformity enhancement is obtained by increasing the number of angular sectors above n = 6.

The analysis of the case of longitudinal slope break helps in understanding how the control strategy acts on the disk orientation. In this case, the difference between the tractor and the ground surface inclination varies in only one direction. Fig. 6 shows the curves of the optimal longitudinal (T) and side (S) tilt signals computed for each disk along the tractor trajectory when 3 angular sectors are considered (n = 3). In this slope situation, the longitudinal tilt signal is the most significant with a maximum value up to 6°. The longitudinal tilt signal is the same for both disks. The side tilt signals are limited to small values less than 1°. This is due to the fact that in the case of a longitudinal slope break, compensation required to keep the particle travelled distances unchanged are more important in the longitudinal direction than in the transverse direction. Furthermore, the side tilt correction curves of both disks are symmetrical because when the right disk is tilted to the left, the left disk is tilted to the right. The longitudinal tilt signal has a bell shape adjusting the orientation of the disk with regards to the inclination of the ground surface. This shape corresponds to the inverted profile of the mean travelled distance of the particles presented in Abbou-ou-cherif et al. (2017a). In the area delimited by x ≤ 5 m, no control is active (T = 0 and S = 0) because the field is a horizontal surface. In the area delimited by x ≥ 30 m, the tractor and the ground surface have the same inclination. Consequently, the longitudinal tilt control is active but with an absolute value less than 1°. The small negative value observed in this situation corresponds to the reduction of the vertical orientation of the particle outlet velocity required to reduce the travelled distances of the particles ejected at the rear of the spreader.

The previous results show that the control strategy based on disks inclination is effective in reducing the application errors under acceptable limits in various slope situations. They provide an overall validation of the control method based on the correction of the particle ballistic flight lengths and of the simplifying assumptions made to solve the optimization problem and reduce the computational complexity. In particular, the study demonstrates that adequate solutions are obtained by using the mean values of particle diameter and ejection parameters; considering only a reduced number of angular sectors; and fitting the 2D cost function through a limited number of points. This is of particular interest to reduce the computation time and thus expect the development of an actual control device working in real-time.

Nevertheless, using the Matlab programming language and a 3.4 GHz quad-core processor, the optimal orientation of the disks could only be derived in less than 0.3 s. This is approximately four times the computation time required for real-time use considering a tractor speed of 12 km/h and a trajectory sampling of 0.25 m. The computation time could be reduced by improving the hardware and using a low level Fig. 4. Digital elevation models (upper row) and corresponding fertilizer application maps simulated without (middle row) or with disk inclination control (lower row) for a 10°longitudinal change of slope (left column), a 10°transverse change of slope (middle column) and an undulating surface with combined slopes set between -6°and +6°(right column). White arrows on the discrete elevation models represent the adjacent passes of the tractor. The rectangular edge drawn on the applications maps delineate the collection areas that are analyzed in Fig. 5. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.) programming language. In order to reduce computation time, the increase of the sampling step of the tractor trajectory could also be considered for further studies.

In this study, the vertical and horizontal mass flow distribution of fertilizer around the disk (and measured experimentally with respect to the disk) was assumed to be unchanging irrespective of the tilt of the disk. Further research is needed to measure potential changes in the mass flow distribution around the disk due to potential changes of the drop point position on the disk or over potential changes of fertilizer motion in the vanes. If the mass flow distribution around the disk changed depending on the disk tilt, then some relationships would be established to describe these changes and introduced in the mechanical model to take into account this behavior.

Conclusion

To sustain fertilizer application uniformity despite disruptions due to the slope and its variations in agricultural field, a new concept of active control of the disk attitude was developed. The control system consisted in offsetting particle travelled distance deviations by changing the inclination of the spinning disk axle. Considering two new degrees of freedom for the disk, the optimal longitudinal and side (i.e. lateral) tilts were derived from a constrained optimization problem. The objective was to minimize the quadratic sum of all the particle travelled distance deviations at each position of the tractor trajectory in the field (i.e. minimizing the variance of the differences between the travelled distances of the particles obtained in the considered topography and the distances that would have been obtained on a horizontal surface). This cost function was simplified by considering a reduced number of angular sectors in the single spread patterns where the flight distance of only one particle weighted by a coefficient can be taken into account. This simplification reduced the computation time and paved the way for the development of real-time control devices. Using a simulation model, three cases of longitudinal, lateral and combined slopes were investigated. Simulations showed that the tilt correcting control reduces effectively the magnitudes of fertilizer application errors, below the acceptable limit of ± 10%, while the application errors could reach 40% when the application maps were simulated without control strategy. Furthermore, the study showed that the choice of a number of angular sectors less than or equal to six is enough to reach relevant disk tilt corrections. Accordingly, simulation results demonstrated the overall performance of the method to achieve adequate uniformity on undulating fields. In practice, the development of control devices could benefit from the development of on-board sensors that are useful to provide model parameters such as the outlet velocity or the mass flow distribution around the disk, in order to reduce preliminary experimental measurements or control these parameters with a closed loop system.

Fig. 1 .

 1 Fig. 1. Kinematic diagram showing the new revolute joints providing two additional degrees of freedom to each disk.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Comparison of two particle travelled distances in the case of a flat field (p 01 and p 02 ) and in the case of a non-flat field (p 1 and p 2 ). (O,x,y,z) is a Cartesian coordinate system.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Mean application rates, with and without disk inclination control, calculated inside the rectangular areas shown on the application maps of Fig. 4 for the 10°l ongitudinal change of slope (a), the 10°transverse change of slope (b) and an undulating surface with combined slopes set between -6°and +6°(c). The mean application rates were computed with various number n of angular sectors used in the optimization problem. (For interpretation of the reference to color in this figure legend, the reader is referred to the web version of this article.)

Acknowledgment

This work was supported by the French National Research Institute of Science and Technology for Environment and Agriculture (Irstea). We would also like to acknowledge the grant of the European Regional Development Fund (FEDER) through the Auvergne Region of France. Finally, we would like to thank the company RAUCH Landmaschinenfabrik GmbH for its support.