Romane Scherrer
email: romane.scherrer@hotmail.fr

Rodrigue Govan

Thomas Quiniou

Thierry Jauffrais

Hugues Lemonnier

Sophie Bonnet

Nazha Selmaoui-Folcher

Automatic Plankton Detection and Classification on Raw Hologram with a Single Deep Learning Architecture

Keywords: Object Detection, Deep Learning, Plankton, Holography

Digital holography is an imaging process that encodes the 3D information of objects into a single intensity image. In recent years, digital holography has been used to detect and count various microscopic objects and has been applied in submersible equipment to monitor the distribution of plankton. To count and classify plankton, conventional methods require a holographic reconstruction to decode the hologram and identify objects. To avoid this iterative and time consuming reconstruction, we demonstrate an object detection based approach that simultaneously performs the detection and classification of all plankton within raw holograms. The proposed method achieved a mean Average Precision (AP@.5) score of 0.981 with 13 plankton species on 512x512 holograms with an inference time of 14 ms per image. Our method can be implemented to detect and count other microscopic objects in raw holograms.

1 Scientific Background The observation and statistics of marine plankton are essential to measure the health of our oceans. In recent years, several submersible equipment [START_REF] Liu | Research progress of digital holography in deep-sea in situ detection[END_REF] (ISIIS, LISST-Holo, eHoloCam) have been deployed as part of large-scale campaigns to acquire in situ images of plankton. Some of these systems use digital holography, a method that enable high resolution images acquisition over a large water column and at high flow rates. Since a hologram encodes the 3D information of all plankton as a single intensity image, a decoding process, called holographic reconstruction, is required to retrieve the sample image from its hologram. Unfortunately, the methods used to process holograms and then count and classify the species are still very time-consuming and manual.

With the multiplication of collected images, various efforts have been made to accelerate and improve the holographic reconstruction, for instance, by adopting a convolutional neural network (CNN) to automatically find the focus [START_REF] Shimobaba | Convolutional Neural Network-Based Regression for Depth Prediction in Digital Holography[END_REF] or to reconstruct a de-focused hologram without performing an auto-focusing or phase recovery routine [START_REF] Rivenson | Phase recovery and holographic image reconstruction using deep learning in neural networks[END_REF][START_REF] Wu | Extended depth-of-field in holographic image reconstruction using deep learning based auto-focusing and phase-recovery[END_REF]. Even though those approaches greatly accelerate the holographic reconstruction, the classification and detection of the objects need to be performed afterwards. As a result, some research works have focused on classifying objects directly on raw holograms, thus avoiding iterative and time-consuming operations required for reconstruction. In [START_REF] Zhang | Automatic classification of marine plankton with digital holography using convolutional neural network[END_REF] a CNN is trained to perform a classification on cropped simulated holograms with 16 different plankton species but a preliminary detection is necessary to determine the regions of interest (ROIs). [START_REF] Trujillo | Automatic detection and counting of phase objects in raw holograms of digital holographic microscopy via deep learning[END_REF] use a sliding window detection method with a CNN that performs a binary classification to detect and count cells. However, all the proposed methods often perform a single task on raw holograms, either a detection or a classification. Object detection methods such as [START_REF] Padilla | A comparative analysis of object detection metrics with a companion open-source toolkit[END_REF] Faster-RCNN, YOLO, SSD or RefineDet, offer an alternative by performing real-time localization and classification of all objects within an image in one pass. By labeling the detected objects within a hologram with rectangular bounding boxes, a detection model could greatly improve the application efficiency of digital holography.

The aim of the paper is to show that the classification and localization of plankton can be simultaneously performed on raw holograms with an object detection model. For that purpose, a dataset of labeled in-line holograms will be simulated with 13 different plankton species. The paper is organized as follows. In Section two, the generation of holographic datasets and the object detection models are described. Section three show the performances of the models. Conclusions and perspectives are given in the last section.

2 Materials and Methods For an in-line configuration (see Fig. 1), the reference and object waves share the same optical axis and an object can be described by a complex transmission function [START_REF] Latychevskaia | Practical algorithms for simulation and reconstruction of digital in-line holograms[END_REF] at a given z plane:

Hologram Formation

t z (x, y) = exp[-a(x, y)] exp[iφ(x, y)] (1)
where a(x, y) describes the absorption of the object and φ(x, y) is the phase distribution.

The transmission function can be used to calculate the wavefront just behind the object

U z + (x, y): U z + (x, y) = t z (x, y)U z -(x, y) (2)
where U z -(x, y) is the incident wave that can be either plane or spherical.

Considering that the object is located at z = 0, the exit wave given by Eq. 2 can be rewritten as U 0 + (x, y) = t 0 (x, u)U 0 -(x, y) and is propagated to the detector/hologram plane which is located at z = z 2 along the optical axis. This propagation is simulated by the angular spectrum method by calculation of the following transformation:

U z 2 (X, Y ) = T F -1 T F (U 0 + (x, y)) × exp 2πiz 2 λ 1 -(λu) 2 -(λv) 2 (3)
where λ is the wavelength and (u, v) are the Fourier domain coordinates. T F -1 and T F denoted the inverse Fourier transform and the direct Fourier transform, respectively.

Note that Eq. 3 is often expressed as

U z 2 (X, Y ) = R(X, Y ) + O(X, Y )
where R and O are the reference and the object waves that interfere at the surface of the recording medium. The recorded hologram at z = z 2 is the intensity calculated by:

H z 2 (X, Y ) = |U z 2 (X, Y )| 2 = U z 2 (X, Y )U * z 2 (X, Y ) (4)
where * denotes the complex conjugate. As a result, a hologram can be simulated once λ, z 2 , U 0 -(x, y) and t 0 (x, y) are known or set.

Transmission Function Simulation

To generate a dataset of labeled holograms for an object detection task, the complex transmission function t 0 (x, y) of several objects in a plane (x, y, z = 0) must be simulated first. In that case, two labeled datasets of plankton images will be used as objects. The first dataset consists of shadow images collected by the In Situ Ichthyoplankton Imaging System (ISIIS) which was the subject of a competition on Kaggle 1 . This open source dataset consists of 121 marine plankton species among which 10 species with a number of images greater than 1000 were selected for our simulations. The second dataset (custom) consists of optical microscopy images of 3 phytoplankton species from New Caledonia (noted P1, P16 and P17). The plankton was imaged with a brightfield microscope at ×10 magnification. The images were thresholded, segmented into (ROIs) and manually labeled. Fig. 2 presents the number of images per species. Note that for each datasets, the ROI segments are labeled and saved as grayscale images where the background has an uniform value equal to 1. As a result, an image I(x, y) can be converted into an absorption a(x, y) = -1 × I(x, y) + 1 so that t 0 (x, y) = exp[-a(x, y)] exp[iφ(x, y)] inside the object support and t 0 (x, y) = 1 where there is no object. In the latter situation, the incident wave remains undisturbed.

To simulate t 0 (x, y) with various objects, the transmission functions of several plankton images can be randomly placed on a N ×N empty (all-ones) image. By doing so, the (x, y)-axis coordinates of the bounding boxes are randomly set. Moreover, the plankton images are already saved as ROIs so that the bounding boxes width and height are the images dimensions. Since the images are classified per species, the labels of a simulated t 0 (x, y) for an object detection task (classes and bounding boxes coordinates) can be completely set.

Hologram Simulation

Before simulating the t 0 (x, y) and the holograms, the plankton images from the two presented datasets were randomly split, per class, in a 80:20 ratio for training and testing, respectively. We have considered that the plankton are pure amplitude objects so that φ(x, y) = 0. The simulation of t 0 (x, y) proceeds as follows. First, for each simulated t 0 (x, y), 13 plankton images (one per species) are randomly selected. The images are then randomly rotated and flipped with 4 possible rotations (0 • ,90 • ,180 • or 270 • ) and 3 possible flips (None, horizontal or vertical). Then, the transmission functions are individually modified so that t plankton (x, y) = exp[-C × a(x, y)] where C is a random constant and C ∈ [0.5, 1]. Next, the 13 transmission functions are randomly placed without overlapping on a 512 × 512 empty image to generate t 0 (x, y). Finally,the hologram H z2 (X, Y ) is simulated with Eq. 3 and Eq. 4. Both the holograms and the t 0 (x, y) are normalized between 0 and 1 and saved.

For the simulations, we have considered λ = 530nm (green), z 2 = 0.8mm and an incident plane wave U 0-(x, y) = 1. 8000 and 2000 holograms were simulated for training and testing, respectively. Fig. 3 presents an example of a simulated and labeled t 0 (x, y) and its corresponding hologram.

Object Detection Models

To perform object detection on raw holograms, we chose two YOLOv5 [START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF] models that were pre-trained on the COCO dataset, namely YOLOv5s2 (the smallest) and YOLOv5x (the largest) with 7.3M and 87.7M parameters, respectively. The models were trained on 8000 holograms during 400 epochs with a batch size of 8 and tested on 2000 holograms. The SGD optimizer was used with an initial learning rate equal to 0.01. To further evaluate the object detection performances on raw holograms, two models were also trained on the transmission functions t 0 (x, u) with the same hyperparameters.

The experiments were carried out on a 2.9 GHz Intel Core i7 PC with 64 GB of RAM and a Nvidia GTX 2060 GPU. The training took 8 hours for the small model and 2 days for the larger one.

Results

In this section, we report the object detection performances with the well-known average precision (AP) metrics [START_REF] Padilla | A comparative analysis of object detection metrics with a companion open-source toolkit[END_REF]. We recall that the AP@.5 and AP@.75 are the average precision computed with an intersection over union (IoU) threshold t = 0.5 and t = 0.75, respectively. The AP@[.5:.95] is reported by computing the mean AP@ with 10 different IoU thresholds [.5:.05:.95]. The mean inference time was evaluated on 2000 holograms and was carried out on a GTX 2060 GPU. It includes FP16 inference, postprocessing and non-maximum suppression. The Tab. 1 summarizes the performances of the object detection tasks performed on the raw holograms and on the transmission functions. For the models trained on the holograms, the AP@.5 are 0.976 and 0.981 for YOLO-v5s and YOLOv5x, respectively. For the models trained on the transmission functions, the AP@.5 are slightly better with 0.985 and 0.993 for YOLOv5s and YOLOv5x, respectively. The AP@[.5:.95] are significantly higher on t 0 (x, y) than on holograms (eg. 0.980 vs. 0.855 for YOLOv5x) but the AP@.75 are still high on holograms (0.928 and 0.955 for YOLOv5s and YOLOv5x, respectively). Those results suggest that the detectors trained on the holograms are efficient for a IoU threshold ≤ .75 but that their performances start to decline at a higher threshold.

Fig. 4 shows the confusion matrix of YOLOv5x at IoU.5 on the test holograms and one of its predictions. One can notice that the diffraction pattern of an object spreads beyond its bounding box. In fact, the further away the object is from the camera, the more this effect will be visible on the hologram. Because of this and the lack of sharp edge, a detector trained on holograms was expected to have difficulty in determining the object boundaries with a high IoU. Although the models trained on t 0 (x, y) perform better than those on holograms, the performances are still promising: the species and location of all the plankton within a 512x512 hologram can accurately be determined with a true positive rate between 91% and 100% for the 13 classes. Moreover, the prediction can be performed in real-time (4 ms for YOLOv5s and 12 ms for for YOLOv5x).

Conclusion

In this paper, we have generated a simulated and labeled holograms dataset from optical and shadow plankton images for an object detection task. 13 different species were considered. Two versions of YOLOv5 are trained on the holograms and on the transmission functions, respectively. Even though the AP of the models trained on the holograms are lower than those trained on the transmission functions, the difference in AP@.5 is less than 1.2%. With a AP@.5 score of 0.981, a YOLOv5x model can perform detection and classification of all plankton within a 512x512 raw holograms in a single pass in 14ms. Although this method was validated with plankton images, it can be implemented to localize, count and identify other microscopic objects in raw holograms.

While the results on simulated holograms are promising, it is often complicated and time consuming to put together a large dataset of real labeled holograms to train a detector. When a small labeled dataset is available, it might be beneficial to pre-train a detector with a large amount of simulated holograms and then use a transfer learning method to fine tune the model on the small dataset. Another approach would be to rely on an intensive data augmentation. Some works in the literature use de-focused backpropagated holograms as inputs of a deep learning model rather than raw holograms. By back-propagated the holograms on several planes near the the correct global focus, the dataset could be significantly enlarged. For future work, we intend to study the object detection performances on a de-focused holograms dataset in which we will address the problem of overlapping objects.

Figure 1 :

 1 Figure 1: In line-holography.

Figure 2 :

 2 Figure 2: Number of images per species. 1 https://www.kaggle.com/c/datasciencebowl/

Figure 3 :

 3 Figure 3: Simulation example. 13 plankton images are used to generate a labeled hologram.

Figure 4 :

 4 Figure 4: Confusion matrix at IoU.5 and model predictions (blue : ground-truth , red : predicted).

Table 1 :

 1 Performances.

	Model	Inputs	AP@.5:.95 AP@.5 AP@.75 Speed
	YOLOv5s	Holograms t 0 (x, y)	0.820 0.967	0.976 0.985	0.928 0.985	4 ms
	YOLOv5x	Holograms t 0 (x, y)	0.855 0.980	0.981 0.993	0.955 0.989	14 ms

https://github.com/ultralytics/yolov5


Holger Döring
email: haw-doering@t-online.de

Occams Razor -No Dark Matter in Galactic Halos

Keywords: dark matter, no dark matter, halo-curves, gravity-forces, galactic core, third Keplerlaw, centripetal force, stellar halo-velocity

Since several years there is the discussion about dark matter in the halos of galaxies. This is a consequence of observation of halo-velocity curves in dependence of distance which don't fulfill the third Kepler-law. Maybe the ansätze are wrong. Shown is a simple physical model in classical Newton-dynamics which fits qualitative the form of observation curves. Additional Parameters can be used to fit the curves quantitative exactly without the hypothesis of dark matter.

1.Introduction:

The analysis of the orbital speeds of stars in spiral galaxies by Vera Rubin since 1960 seem to show the problem: The orbital speed of the stars in the halos of galaxies would have to be much lower than it actually is with increasing distance from the galaxy's center after the classical Newtonian model resp. the third law of Kepler. So since several years the existence of a form of "dark matter" is discussed in the scientific literature with some very fantastic explanations and speculations about the reason of this observation like is usual in modern theoretical physics, astrophysics and astronomy [1.],[2.],[3.],[4.],[12.] A few mathematical models are created to clear this problem, some of them modified Newton Dynamics (MOND), Yukawa-like exponential developments or added discussed gravity-effects caused by GRT in large galactic scales ( gravity-caused time-dilation at the edge of the galaxies) etc.[5],[6],[7], [8]. Also some Scalar-vector-tensor theories are invented and mentioned (and permutated). But if the principle of Occams Razor in physics is taken seriously there should be a form of trying a short and easy explanation of this phenomenon without going too far into fantastic elaborated models. Sometimes the discussion of this problem looks like the aether-discours in nineteenth century before Poincaré, Lorentz and Einstein ended this by solving it via formulation of SRT. (But on the other hand Ernst Mach denied the existence of atoms, because at his time the experimental data wasn't sufficient, there has to be caution in argumentation ) and "dark matter" seems like a spooky ghost of rumour, which was never seen but all speak about it like aether in the 19 th century.

The summed rotation curve over radius of a galaxy adds from galactic gas, from the galactic disc and the halo-curve [9.],[10.], [11.]. Only the halo observations causes problems with theory in a galaxy so this discussion here is only about the halo-problem. Dark matter between galaxies in greater clouds (observed by Zwicky and later) may be exist and therefore about this form of DM there are no statements in this paper here. Other possible forms of DM, which are used in explaining phenomena at gravitational lenses are not discussed either. Both forms may exist or not. This paper here is only about galactic states of stellar rotation curves.

Example given for the halo-problem in the following picture:

2.The common thinking:

The common explanations are: setting centripetal force equal to gravitational force.

F Z =F G
This leads directly to Keplers third law:

v = √ G⋅M r .
The galactic disc fulfills this law approximatively but the halo doesn't.(For the planets in our solarsystem its very exactly fulfilled). Observed is a form of v ∼r or v ∼r 1/ 2 instead of the expected

Graph 1:In the picture above shown is the galactic disk which can be described via Keplers 3 rd law and "dark matter"means the galactic halo, which curve causes the problem of difference between observation and explanation.

v ∼r -1/ 2 . The conclusions are: therefore there has to be a form of "dark matter" in the halos of the galaxies , which causes this effect. Shown is now, that this assumption is unnecessary with the principle of Occams Razor to take the nearest explanation instead of going fantastic like it is modern. Therefore there are now calculations made in the following text with normal Newton Law to fulfill the requirements:

3.Calculation:

2.1. The gravitational Force-model:

Forces are now positioned like ⃗ F G ,C -⃗ F G , H = ⃗ F Z , H
with the following force-diagrams as reason why.

Graph 2: The vector -diagram of effective forces working at galactic core and halo.

Explanation: F(G,C) is the gravitaional force at the galactic core. F(G,H) is the gravitational force in the halo and F(Z,H) is the centripetal force of the halo.

Therefore there is:

M⋅m⋅G R 2 - M⋅m⋅G (R+r ) 2 = m⋅v ² R+ r (1.)
This leads directly to:

Graph 3: Shown is a simple model of galactic core with vectors of F(G,C) and galactic halo with F(Z,H) and F(G,H)

.

v ( r)= √ M⋅G⋅ ( R+r R ² - 1 R+r ) (2.)
Here R is core-radius, measured from galactic-center and r is halo-radius, measured from galactic core. R+r is whole halo-radius from galactic center. Answer: the model-curve looks qualitatively like the observed halo-curve. But in the limit for r →∞ there is no constant value of v(r) but it diverges weak to infinity. So the model has probably to be rejected because it probably doesn' fit with edges of observation curves.But it is a near approximation for r not to great.

Result

The energy-model for kinetic and potential energies:

From forces to energy there is:

E kin = | ∫ 0 R +r F Z ,H dr | = | ∫ R R +r F G , H dr | = E Pot (3.)
This leads directly to:

v H = √ 2⋅M⋅G⋅ ( 1 R - 1 R+r ) -v C (4.)
If the left side of F Z is integrated from R to R+r , then v c vanishes to zero. This ansatz above is the simplest way of calculating this problem but there are some three others more complicated. All three formulas lead to the same result excepting one which added a numerical factor of 2 in R .If the curve is plotted in v (r) , there will be the following graph:

Result: the velocities convergate strong and go straight into a constant value for lim r →+ ∞ v (r ) .

This constant limit for the home galaxy would be v H (r )= √ 2⋅M⋅G R ≈557,7 km s .

If the sun would be in the halo, the energy-model produces a hypothetical value of v (8,5 kp)≈443 km s with same constants as above.

3.Conclusion:

The force-model has probably to be rejected because the velocities are too strong in increasing and have no constant limit like is observed in reality for halo star-systems. The energy-model may be a good explanation for star moving behaviour in the halos of galaxies because it fits qualitatively with observation-curves. It is obvious, that exact measurements from some galaxie-halos are necessary to test the model and to take it seriously or to reject it.

4.Summary:

Two models are put up to explain the moving of halo-stars in galaxies without dark matter. One of the model has probably to be rejected, the force-model, but the other model fulfills at least the quality of the demanded observation laws. Further exact details from measuring values of some Graph 5: Qualitative graph of v (r) in km/s from energy conservation-conditions with M=0.9*200*10 sunmasses for the home-galaxy in the core and distance r in kpc. ⁹ galaxies have to be made to prove the value of this energy-model and to fit it in with observation or to reject it either.

4.References:

  : there is no halo mass m in this equation (2.). All movements are determinded by core mass M . If the curve of equation (2.) is plotted , there is the result of picture 4: Example given: Numerical data for this force-model: If the sun, which is in the disc ( v S , D =220 km s R=5 kp ; r=r S =8,5 kp ; M =0,9 * 2 * 10¹¹ * 2 * 10³⁰ kgR -core radius from galactic center, r -suns radius from galactic center M -90% from galactic mass in sun masses (kg) in core supposed (classical case).

Graph 4 :

 4 this force-curve is plotted for M=0.9*200*10 sunmasses in the core as function of v(r) ⁹ over r. Distance r in kpc. Velocity v in km/s.


Philippe Grangier
email: philippe.grangier@elementai.com

Michel Gendreau
email: michel.gendreau@cirrelt.net

Fabien Lehuédé
email: fabien.lehuede@imt-atlantique.fr

Louis-Martin M Rousseau
email: louis-martin.rousseau@cirrelt.net

The vehicle routing problem with cross-docking and resource constraints

In this paper, we propose an extension of the vehicle routing problem with cross-docking that takes into account resource constraints at the cross-dock. These constraints limit the number of docks that can be used simultaneously. To solve this new problem, we adapt a recently proposed matheuristic based on large neighborhood search. In particular, we focus on the feasibility tests for insertions and compare heuristics and constraint programming strategies. Finally, computational experiments on instances adapted from the vehicle routing problem with cross-docking are reported. They give insights on the impact of a limited cross-dock capacity on the routing cost.

Introduction

In logistics, cross-docking is a distribution strategy in which goods are brought from suppliers to an intermediate transshipment point, the so-called cross-dock, where they may be transferred (without storing) to another vehicle for delivery. Compared to traditional distribution systems, cross-docking can help reducing delivery costs and delivery lead time, that is why it is used by many companies from different sectors: less-than-truckload (LTL), retail or automotive for example [START_REF] Van Belle | Cross-docking: state of the art[END_REF]. In the vehicle routing literature, the associated routing problem is called the Vehicle Routing Problem with Cross-Docking (VRPCD) [START_REF] Wen | Vehicle routing with crossdocking[END_REF]. It is a variant of the Pickup and Delivery Problem with Transfers [START_REF] Cortés | The pickup and delivery problem with transfers: Formulation and a branch-and-cut solution method[END_REF] with one compulsory transfer point: vehicles start by collecting items, then return to the cross-dock where they unload/reload some items and eventually visit delivery locations. A few authors have proposed methods to solve this problem or variants of it, but to our knowledge, in most models there is no limit on the processing capacities of the cross-dock: as soon as a truck arrives, it immediately undergoes consolidation operations. However in practice, this may not be the case because of limited equipment or workforce, thus trucks may have to wait before being unloaded. In fact, a wide range of the cross-docking literature is dedicated to the scheduling of operations at the cross-dock taking into account the cross-dock capacity. It has recently been pointed out [START_REF] Buijs | Synchronization in cross-docking networks: A research classification and framework[END_REF][START_REF] Ladier | Cross-docking operations: Current research versus industry practice[END_REF], that there is a need to consider the synchronization of local and network-wide cross-docking operations, in particular to take into account resource capacity at the cross-dock. To that end, in this paper, we introduce a new variant of the VRPCD in which the number of vehicles that can simultaneously be processed at the cross-dock is limited. We call it the Vehicle Routing Problem with Cross-Docking and Resource Constraints (VRPCD-RC). The dock resource constraint is a resource synchronization constraint as defined by [START_REF] Drexl | Synchronization in vehicle routing-A survey of VRPs with multiple synchronization constraints[END_REF] as vehicles compete to access a scarce resource: the processing capacity of the cross-dock. Very often resource synchronization constraints imply a difficult scheduling problem which is embedded within the vehicle routing problem. VRPCD-RC is no exception and a major contribution of this paper is on the integration of the scheduling problem associated with the dock resource constraints within a recently proposed large neighborhood search (LNS) based matheuristic [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD.

The remainder of this paper is organized as follows. A literature review is presented in Section 2, while the problem is defined in Section 3. In Section 4, we present the proposed method. Computational results are presented in Section 5, and a conclusion is given in Section 6.

Literature review

In this section we review the literature on two related vehicle routing problems: the vehicle routing problem with cross-docking and vehicle routing problems with resource synchronization.

The vehicle routing problem with cross-docking

Many cross-docking problems have received attention from operations research and management science researchers. For instance: location, assignment of trucks to doors, inner flow optimization or routing [START_REF] Van Belle | Cross-docking: state of the art[END_REF]. In particular, the vehicle routing problem with cross-docking consists in designing routes to pick up and deliver a set of transportation requests at minimal cost using a single cross-dock. It was introduced by [START_REF] Lee | Vehicle routing scheduling for cross-docking in the supply chain[END_REF] in a variant which imposes trucks to arrive at the exact same time at the cross-dock. [START_REF] Wen | Vehicle routing with crossdocking[END_REF] relaxed this constraint only imposing precedence constraints based on the consolidation decisions and they also added time windows. This is the most studied variant, and it is the one we will refer to as the vehicle routing with cross-docking (VRPCD). Several heuristics have been proposed to solve it: based on tabu-search [START_REF] Wen | Vehicle routing with crossdocking[END_REF][START_REF] Tarantilis | Adaptive multi-restart tabu search algorithm for the vehicle routing problem with cross-docking[END_REF][START_REF] Nikolopoulou | Moving products between location pairs: Cross-Docking versus Direct-Shipping[END_REF][START_REF] Nikolopoulou | Adaptive memory programming for the many-to-many vehicle routing problem with crossdocking[END_REF], iterated local search (Morais et al., 2014) and LNS [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF]. Many variants have been studied, with optional cross-dock returns [START_REF] Petersen | The pickup and delivery problem with cross-docking opportunity[END_REF][START_REF] Santos | The pickup and delivery problem with cross-docking[END_REF][START_REF] Nikolopoulou | Moving products between location pairs: Cross-Docking versus Direct-Shipping[END_REF], taking into account the cost of handling (Santos et al., 2011a,b;[START_REF] Enderer | Integrating dock-door assignment and vehicle routing with cross-docking[END_REF], with multiple cross-docks [START_REF] Ahmadizar | Two-level vehicle routing with crossdocking in a three-echelon supply chain: A genetic algorithm approach[END_REF][START_REF] Maknoon | Vehicle routing with cross-dock selection[END_REF][START_REF] Tation Kroep | The Two-period Vehicle Routing Problem with Crossdocking and Transfers[END_REF], with partly subcontracted or left to suppliers routes [START_REF] Yu | Open vehicle routing problem with crossdocking[END_REF][START_REF] Enderer | Integrating dock-door assignment and vehicle routing with cross-docking[END_REF]. [START_REF] Dondo | A monolithic approach to vehicle routing and operations scheduling of a cross-dock system with multiple dock doors[END_REF] considered a case where the number of doors at the cross-dock is fixed and smaller than the number of trucks. In particular each door is modeled individually: a time matrix models the time spent by a truck for moving from an inbound door a to an outbound door b. They solved two randomly generated instances with up to 70 requests with a mathematical model combined with a sweep heuristic.

From a general perspective, the VRPCD can be viewed as a special case of the pickup and delivery problem with transfers with only one compulsory transfer point [START_REF] Cortés | The pickup and delivery problem with transfers: Formulation and a branch-and-cut solution method[END_REF]. [START_REF] Guastaroba | Intermediate facilities in freight transportation planning : a survey[END_REF] released a survey on intermediate facilities in freight transportation. For a general overview of cross-docking and crossdock related problems we refer the reader to [START_REF] Boysen | Cross dock scheduling: Classification, literature review and research agenda[END_REF]; [START_REF] Agustina | A review: Mathematical models for cross docking planning[END_REF][START_REF] Van Belle | Cross-docking: state of the art[END_REF][START_REF] Nassief | A comparison of formulations and relaxations for cross-dock door assignment problems[END_REF] for a recent article on the dock-door assignment problem.

Resource synchronization

The expression resource synchronization appears in [START_REF] Drexl | Synchronization in vehicle routing-A survey of VRPs with multiple synchronization constraints[END_REF] as a way to model the following constraint: 'The total consumption of a specified resource by all vehicles must be less than or equal to a specified limit.' Of course, this resource has to be scarce to be constraining, as vehicles compete to access it. Such resource constraints arise in many different vehicle routing problems usually when a special infrastructure or equipment is required: a docking station or parking space in airport cargo system [START_REF] Ebben | Dynamic transport scheduling under multiple resource constraints[END_REF], a berth in maritime transportation [START_REF] Gronhaug | A Branch-and-Price Method for a Liquefied Natural Gas Inventory Routing Problem[END_REF], a forest loader in forestry [START_REF] Hachemi | A hybrid constraint programming approach to the log-truck scheduling problem[END_REF], a pump in ready mix-concrete delivery [START_REF] Schmid | A Hybrid Solution Approach for Ready-Mixed Concrete Delivery[END_REF], an asphalt paver in public works [START_REF] Grimault | A two-phase heuristic for full truckload routing and scheduling with split delivery and resource synchronization in public works[END_REF][START_REF] Grimault | An adaptive large neighborhood search for the full truckload pickup and delivery problem with resource synchronization[END_REF], intermediate charging stations for electric vehicle [START_REF] Enderer | Integrating dock-door assignment and vehicle routing with cross-docking[END_REF]. Limited storage [START_REF] Ebben | Dynamic transport scheduling under multiple resource constraints[END_REF] or processing capacities [START_REF] Hempsch | Vehicle routing problems with inter-tour resource constraints[END_REF] can also account for resource synchronization constraints. [START_REF] Hempsch | Vehicle routing problems with inter-tour resource constraints[END_REF] also mention a situation where only a fixed number of vehicles (smaller than the total number) can perform long routes. When the constrained resource travels to deliver products or services, it corresponds to the definition of resource constrained routing given in the recent survey by [START_REF] Paraskevopoulos | Resource constrained routing and scheduling: Review and research prospects[END_REF].

Many approaches have been applied to deal with these resource constraints in vehicle routing problems. [START_REF] Ebben | Dynamic transport scheduling under multiple resource constraints[END_REF] sequentially inserted requests and check resource constraints in a predefined order. El [START_REF] Hachemi | A hybrid constraint programming approach to the log-truck scheduling problem[END_REF] use a dedicated constraint programming model, later combined with a greedy scheduling heuristic in El [START_REF] Hachemi | A heuristic to solve the synchronized log-truck scheduling problem[END_REF], to ensure that resource constraints are satisfied during their entire solving process. In ready-mix concrete routing problems as well as in public works routing problem, orders are larger than truck capacity, as such they have to be split into several delivery operations that should not overlap. Resource synchronization constraints arise at pickup sites or at delivery sites or at both. [START_REF] Asbach | Analysis, modeling and solution of the concrete delivery problem[END_REF], [START_REF] Schmid | A Hybrid Solution Approach for Ready-Mixed Concrete Delivery[END_REF], [START_REF] Schmid | Hybridization of very large neighborhood search for ready-mixed concrete delivery problems[END_REF] and [START_REF] Grimault | An adaptive large neighborhood search for the full truckload pickup and delivery problem with resource synchronization[END_REF] imposed precedence constraints on the sequencing of operations to handle precedence constraints. [START_REF] Gronhaug | A Branch-and-Price Method for a Liquefied Natural Gas Inventory Routing Problem[END_REF] relied on a time-discretized formulation in which resource synchronization constraints are easily expressed. [START_REF] Hempsch | Vehicle routing problems with inter-tour resource constraints[END_REF] proposed a generic modeling for inter-routes constraints via the use of Resource Extension Functions (REF). In particular they focused on the use of REF as efficient feasibility tests in local search based algorithms when the solution is represented by a giant tour. On the other hand, Froger et al. ( 2017) relax charging station capacity constraints when generating routes via an iterated local search for an electric vehicles routing problem. Only when assembling routes into a solution do they enforce these synchronization constraints in a Benders' like approach. They call this approach a route-first assemble-second. Their later part shares some similarity with the approach that is proposed in this paper, however in our method we enforce synchronization constraints even during the route generation phase. Generally speaking, these type of approaches are column generation based matheuristics in the framework of [START_REF] Archetti | A survey on matheuristics for routing problem[END_REF].

From this literature review, it is clear that most of the resource synchronization constraints are in fact complex scheduling problem integrated into a vehicle routing problem. However the precise definition of the scheduling problem depends largely on the vehicle routing problem at stake.

The vehicle routing problem with cross-docking: model and resource synchronization constraint

In the VRPCD-RC, we consider a set of transportation requests, each having different origins and destinations. Each vehicle starts from the depot, performs a pickup leg which ends at the cross-dock, then it performs a delivery leg and returns to the depot. At the cross-dock, some transportation requests are exchanged between vehicles, creating precedence constraints between pickup and deliveries legs of different vehicles. We consider the number of doors that can be processed simultaneously a the cross-dock is limited. This new constraint may introduce some additional waiting times for vehicles at the cross-dock (compared to the VRPCD). We consider two cross-dock configurations:

-In the shared configuration, the total number of vehicles that can be processed simultaneously is limited to a number S. -In the separated configuration, the unloading capacity is separated from the loading capacity into inbound and outbound doors.

This section presents the vehicle routing problem with cross-docking, and in particular the cross-dock model, as defined by [START_REF] Wen | Vehicle routing with crossdocking[END_REF]. It also presents in details the two considered cross-dock models.

The vehicle routing problem with cross-docking

In the VRPCD, we consider a cross-dock c, a set of requests R, and a homogeneous fleet of vehicles V , each of capacity Q and based at o. Each request r ∈ R has to be picked up at its pickup location p r within its pickup time window [e pr , l pr ], and delivered at its delivery location d r within its delivery time window [e dr , l dr ]. In the case of early arrival, a vehicle is allowed to wait, but late arrivals are forbidden. We denote by P the set of pickup locations and by D the set of delivery locations. Each vehicle starts at o, goes to several pickup locations, and then arrives at the crossdock where it unloads/reloads some requests. It then visits several delivery locations and eventually returns to o. Note that a vehicle must visit the cross-dock even if it does not unload nor reload there. The sequence of operations at the cross-dock is described in Section 3.2. The base of vehicles o and the cross-dock c may be the same location. The VRPCD is defined on a directed graph G = (V, A), with G = {o}∪P ∪{c}∪D and

A = {(o, p)|p ∈ P } ∪ P × P ∪ {(p, c)|p ∈ P } ∪ D × D ∪ {(d, o)|d ∈ D} ∪ {(o, c), (c, o)}.
With each arc (i, j) ∈ A is associated a travel time t i,j and a travel cost c i,j . Solving the VRPCD involves finding |V | routes, and a schedule for each route, such that the capacity and time-related constraints are satisfied, at minimal routing cost. An arc-based mathematical formulation can be found in [START_REF] Wen | Vehicle routing with crossdocking[END_REF].

Precedence constraints at the cross-dock

Following [START_REF] Wen | Vehicle routing with crossdocking[END_REF], if a vehicle k has to unload a set of items R - k and reload a set R + k at the cross-dock, the time spent at the cross-dock can be divided into four periods:

-Preparation for unloading. The duration δ u of this period is fixed.

-Unloading. The duration of this period depends on the quantity of items to unload.

For vehicle k the duration is ( i∈R - k q i )/s u , where s u is the unloading speed in quantity per time unit. All unloaded items become available for reloading at the end of this period.

-Preparation for reloading. The duration δ r of this period is fixed.

-Reloading of requests. The duration of this period depends on the quantity of items to reload. For vehicle k the duration is ( i∈R + k q i )/s r , where s r is the reloading speed in quantity per time unit. All the items for loading must have been unloaded before the beginning of the reloading operation (preemption is not allowed).

Items that are not transferred at the cross-dock remain in the vehicle. Thus, if a vehicle does not unload or reload it need not spend any time at the cross-dock and can leave immediately.

Resource constraints models at the cross-dock

According to [START_REF] Van Belle | Cross-docking: state of the art[END_REF], the most common service mode of a cross-dock is called exclusive. In an exclusive mode, a dock door is either exclusively dedicated to unloading (inbound operations) or reloading (outbound operations). Such assignment is a decision taken at a strategical or tactical level that cannot be modified in the VRPCD, which is an operational problem. In practice most cross-docks are I-shaped [START_REF] Van Belle | Cross-docking: state of the art[END_REF], with inbound doors on one side and outbound doors on the other. The flow of items is thus uni-directional, which is typically easier to manage. This is a common situation but it is not mandatory.

Cross-docks usually have many doors (typically ranging from 40 to 150 according to [START_REF] Van Belle | Cross-docking: state of the art[END_REF]). However, as mentioned in [START_REF] Hempsch | Vehicle routing problems with inter-tour resource constraints[END_REF] or [START_REF] Li | Crossdocking -JIT scheduling with time windows[END_REF], processing capacities may actually be lower than the number of doors. This comes from a limited workforce or special equipment to move the items within the cross-dock. We denote by A I (resp. A O ) the number of inbound doors (resp. outbound) that are being processed simultaneously. Because of the previous two considerations regarding operations at the cross-dock, we will consider two exclusive cross-dock models that integrate resource constraints:

a case in which the total number of doors (inbound plus outbound) that can be processed simultaneously is limited to a number S, ie A I + A O ≤ S. This is a simple way to model resource constraints due to limited workforce. We refer to this case as shared. a case in which the number of inbound doors (resp. outbound doors) that can be processed simultaneously is limited to a number

I (resp O), ie A I ≤ I (resp. A O ≤ O)
. This is a simple way to model a resource constraints due to special equipment on each side of the cross-dock. We refer to this case as separated.

Provided that at least two dock doors can be processed simultaneously, the scheduling problems at the cross-dock are NP-Hard as the scheduling problem P 2||C max is included in them. In the rest, we will simply use dock to refer to the capacity in the number of docks processed simultaneously.

We call the VRPCD with these dock resource constraints, vehicle routing problem with cross-docking and resource constraints (VRPCD-RC). Figure 1 illustrates the sequence of consolidation operations for a vehicle that unloads and reloads items at the cross-dock in the VRPCD-RC. When the vehicle arrives at the cross-dock, it can immediately be prepared for unloading, but it has to wait before being actually unloaded (W U ) because of a lack of available resources (note that such waiting time does not exist in the VRPCD). Once the unloading operation is done, the vehicle can proceed and move to an outbound door. Again the preparation can be performed immediately, but the truck may wait before being actually reloaded (W R ). This waiting time can have two origins: first, not all items are available when it is ready (such situation can also arise in the VRPCD), second, there maybe a lack of available resources (which cannot occur in the VRPCD). In this section, we present the matheuristic we propose for the VRPCD-RC: Large Neighborhood Search+Set Partitioning and Scheduling (LNS+SPS). It is derived from the method of [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD. We recall the structure of this previous method in Section 4.1, while the overall LNS+SPS method is summarized in Section 4.2. In Section 4.3 we present the destruction and repair operators used in LNS. The resource constraints require special feasibility test for insertions in the repair methods that are detailed in Section 4.4. The set partitioning component called Set Partitioning and Scheduling problem (SPS) is presented in Section 4.5.

4.1 Matheuristic of [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD In a previous paper [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], we proposed a matheuristic to solve the VRPCD, that is still, to our knowledge, one of the best method for this problem. This method is a LNS with periodic calls to a set partitioning based problem that assembles parts of previously obtained solutions. It is called LNS + SPM (Large Neighborhood Search + Set Partitioning and Matching) LNS [START_REF] Shaw | Using constraint programming and local search methods to solve vehicle routing problems[END_REF] iteratively destroys (removes several requests from) and repairs (reinserts requests into) the current solution using heuristics. In [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], most of these heuristics are adapted to the VRPCD from [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF]. For each request to reinsert, potential transfers lead to quadratic repair possibilities in the number of vehicles (as opposed to linear in the VRP). We only considered transferring a request among vehicles that drive near its pickup or delivery location as a way to fight this growth. This help speeding up the whole method without loss in quality. On top of that, by reusing the technique of [START_REF] Masson | Efficient feasibility testing for request insertion in the pickup and delivery problem with transfers[END_REF] time feasibility of reinsertion could be checked in constant time.

Each time a solution is obtained by applying the destroy+repair process of the LNS, for each route its pickup and delivery parts (called legs) are stored in a memory. The SPM component aims at assembling all the legs in memory into a potentially better solution. Due to possible transfers at the cross dock, there exists some precedence constraints between pickup and delivery legs. This makes the SPM a set partitioning problem with additional precedence constraints. It is worth noting that coupling a set partitioning component and a meta-heuristic in a route-first assemble second fashion has proved being successful for many VRP (see e.g. [START_REF] Rochat | Probability diversification and intensification in local search for vehicle routing[END_REF]; [START_REF] Archetti | An optimization-based heuristic for the split delivery vehicle routing problem[END_REF]; [START_REF] Subramanian | A hybrid algorithm for a class of vehicle routing problems[END_REF]; [START_REF] Pillac | A parallel matheuristic for the technician routing and scheduling problem[END_REF]; [START_REF] Villegas | A matheuristic for the truck and trailer routing problem[END_REF]; [START_REF] Mendoza | A hybrid metaheuristic for the vehicle routing problem with stochastic demand and duration constraint[END_REF][START_REF] Enderer | Integrating dock-door assignment and vehicle routing with cross-docking[END_REF]; [START_REF] Tellez | The fleet size and mix dial-a-ride problem with reconfigurable vehicle capacity[END_REF]).

This algorithmic approach has been retained for the VRPCD-RC, in what follows we will detail each part of the algorithm, specifying if it has been reused from the VRPCD or what changes were made for the VRPCD-RC and why.

Structure of the proposed method

Algorithm 1 presents a sketch of the proposed method: Large Neighborhood Search + Set Partionning and Scheduling. The main component is LNS (l.5-16), which is enhanced by the occasional resolution of an SPS (l.20). When the LNS finds a new solution, both routes (l. 18) and pickup and delivery part of each route (l. 17) are added to memory components called pools. The SPS (Algorithm 2) aims to find the best possible solution from the route in the pool of routes after applying some dominance-based replacement using the pool of legs. This is detailed in Section 4.5. Hereafter, the destruction and repair methods used in the LNS are presented.

Destruction operators

When partially destroying a solution, a destruction method M -and a number Φ of requests to remove are selected. Unless stated otherwise, this method is reused until Φ is reached. The destruction operators are taken from [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF].

Random removal: a request is removed at random.

Worst removal: a request with a high removal gain is removed, where removal gain is the difference in the cost of the solution with and without the request. The requests are sorted in non-increasing order of removal gain and put in a list N . The request to remove is selected in a randomized fashion as in [START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]: given a parameter p, a random number y between 0 and 1 is drawn. Then the request in position y p × |N | is removed.

Historical node-pair removal: each arc (u, v) ∈ G is associated with the cost of the cheapest solution it appears in (initially this cost is set to infinity). The request that is served using the arcs with the highest associated costs is removed using a randomized selection similar to that of worst removal.

Related removals: these methods aim to remove related requests. Let the relatedness of requests i and j be R(i, j). Two distinct relatedness measures are used: distance and time. The distance measure between two requests is the sum of the distance between their pickup points and the distance between their delivery points. The time measure is the sum of the absolute difference between their start of service at their pickup points and the absolute gap between their start of service at their delivery point.

In both cases a small R(i, j) indicates a high relatedness. A randomized selection, similar to worst removal (albeit with a non decreasing ordering), is performed.

Transfer removal: for each pair of routes (v i , v j ), with v i = v j the number of requests transferred from v i to v j is computed. Then a roulette-wheel selection is applied on the pairs of routes (the score of a pair being the number of requests transferred), and the requests that are transferred between the routes in the selected pair are removed.

If there are less transferred requests than the target number Φ to remove, the rest of the removals is performed with random removal.

Repair operators

In LNS, unplanned requests are stored in a request bank. In the following, we describe the operators that are used to reinsert them in an incomplete solution.

Best insertion: among all requests r in the request bank, the request with the cheapest insertion (considering insertions with and without transfer) is selected and inserted at its best position.

Regret Insertion: for each request r in the request bank and for each pair of vehicles (pickup vehicle, delivery vehicle), the cost of its cheapest feasible insertion (if any) is computed. Note that the pickup and delivery vehicles may be the same (in the case of insertion without transfer). Then, with these insertion options, the k-regret value of r is defined as

c k r = k i=1 (f i -f 1 )
, where f 1 is the cost of the cheapest insertion, f 2 is the cost of the second-cheapest insertion and so on, and k is a parameter. The cheapest insertion of the request with the highest regret value is performed.

For repair operators, two facts should be taken into account: first, because of resource constraints, only a limited number of requests will be transferred (some requests that were transferred in the VRPCD may no longer be transferred in the VRPCD-RC). Thus, there is an incentive in creating routes without transfer. Second, as detailed in Section 4.4.2, feasibility tests can no longer be performed in constant time. Thus, repair methods that check many insertions, such as k-regret with a high value of k, should be avoided. As such we use: best insertion, 2-regret, best insertion without transfer and 2-regret without transfer as repair methods. In the variants without transfer, only insertions without transfer are considered. This is a difference with [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] where all repairs methods were considering transfers and where 3 and 4-regret insertions were computationally affordable.

The initial solution is obtained by applying a 2-regret without transfer.

Integration of dock resource constraints in LNS

In the repair methods of LNS, we need to ensure that the insertions we consider are feasible both with respect to capacity and time-related constraints (time windows and dock resource). Capacity constraints can easily be checked in constant time. In this section we focus on how to handle dock resource constraints. To that end, we start by presenting the scheduling model associated with dock resource constraints, then the methods we propose to solve it and eventually the general structure of feasibility tests we use for maximal efficiency.

Scheduling problems associated with dock resource constraints

For each route k ∈ K, let a k be its earliest feasible arrival time at the cross-dock and b k its latest feasible departure time from the cross-dock. Given an insertion that would insert request r in route k 1 for pickup and route k 2 for delivery, we can compute the new earliest feasible arrival time at the cross-dock of k 1 : a 1 , and the new latest feasible departure time of k 2 : b 2 (provided that no time windows violation occurs in the pickup leg of k 1 and in the delivery leg of k 2 , in which case we could immediately reject the insertion). Thus, dock resource constraints can be seen as a satisfaction scheduling problem at the cross-dock.

For each route k ∈ K, let T k be the set of routes that deliver at least one request picked up by k; let t - k and t + k be its associated unloading and reloading operations respectively; and let R - k and R + k be the sets of requests being unloaded and reloaded respectively. Let isActive be an indicator function such that isActive(o, h) is equal to 1 if and only if task o is being performed at instant h. Let H be the time horizon of the problem. In the separated case the associated scheduling problem is:

t + k .Start ≥ t - k .End ∀k ∈ K, k ∈ T k (1) t + k .Start ≥ t - k .End + δ r ∀k ∈ K; s.t. R + k = ∅ (2) t + k .Start ≥ t - k .End ∀k ∈ K; s.t. R + k = ∅ (3) t - k .Start ≥ a k + δ u ∀k ∈ K; s.t. R - k = ∅ (4) t - k .Start ≥ a k ∀k ∈ K; s.t. R - k = ∅ (5) t + k .End ≤ b k ∀k ∈ K (6) k∈K isActive(t - k , h) ≤ I ∀h ∈ [0, H] (7) k∈K isActive(t + k , h) ≤ O ∀h ∈ [0, H] (8)
In the shared case the associated scheduling problem is:

(1 -6) k∈K isActive(t - k , h) + isActive(t + k , h) ≤ S h ∈ [0, H] (9)
Constraints (1) ensure that all reloading operations that depend on a route k start no earlier than the end of the unloading task associated with k. Constraints [START_REF]and best solution found in the shared cross-dock configuration case[END_REF] and (3) ensure the delay (preparation) between the unloading and reloading task of a route k respects the model presented in Section 3.2. Constraints (4) and ( 5) ensure that for each route, its corresponding unloading operation cannot start before the earliest feasible arrival time at the cross-dock. Constraints (6) ensure that for each route, its corresponding reloading operation is done by its latest feasible departure time. Constraints (7 and 8) model the separated case while constraints (9) models the shared case.

Constraints on start and end of some events as well as the isActive function are easily modelled respectively with the concept of IntervalVar and CumulExpr from the OPL modelling language [START_REF] Van Hentenryck | The OPL Optimization Programming Language[END_REF], which can be used by IBM CP Optimizer (IBM Corporation, 2014). In particular, we recall that (from IBM Corporation (2014)):

'An interval variable represents an interval of time during which a task happen, and whose position in time is an unknown of the scheduling problem. An interval is characterized by a start value, an end value and a size. (...) An interval variable can be optional, that is, one can decide not to consider [it] in the solution schedule.'

Proposed methods for the dock resource constraints

To solve the satisfaction scheduling problems of Section 4.4.1 we propose two methods:

(1) using a third party CP solver or (2) using scheduling heuristics. When repeatedly calling a third party solver we cannot neglect its overhead (model building, memory allocation of the solver, ...) potentially leading to very high run-times. Scheduling heuristics are potentially faster, we could thus perform more LNS iterations within the same time budget, but they are likely to report false negatives (stating that an insertion is unfeasible although it is actually feasible). This is a common run-time versus solution-quality trade-off situation.

The scheduling heuristics we use are list heuristics: among a list of available tasks (a task is said to be available for scheduling if all its predecessors have already been scheduled) we select one task according to a given criterion and we try to schedule it. If a feasible schedule is found, we update the resource constraints accordingly, then we update the list of available tasks to schedule. We repeat this procedure until no more tasks have to be scheduled. If at one point we cannot schedule a task, we declare this insertion unfeasible according to this scheduling heuristic. The four different selection criteria we use are listed hereafter.

First Come First Served (FCFS): the task with the earliest release date in the list of available tasks is selected.

Earliest Due Date (EDD): the task with the earliest due date in the list of available tasks is selected.

Most Successors First (MSF): the task with the largest number of successors in the list of available tasks is selected. By definition reloading tasks do not have successors, we break ties with the EDD rule.

Shortest Processing Time First (SPTF): the task with the shortest processing time in the list of available tasks is selected.

Checking the feasibility of an insertion in the VRPCD-RC

A necessary condition for an insertion to be feasible in the VRPCD-RC is that it is feasible in the VRPCD. As mentioned in Section 4.1, feasibility tests in the VRPCD can be performed in a constant time. On the other hand, feasibility test with respect to dock resource constraints in the VRPCD-RC cannot be done in a constant time. As a result, we test the feasibility of an insertion as shown on Fig. 2: we start by checking if it is feasible for the VRPCD, if the insertion passes this test, we test it with respect to dock resource constraints using either a CP solver or scheduling heuristics. Given a set of routes K, set partitioning and scheduling (SPS) aims to select a subset K of K minimizing cost such that (1) each request is picked up and delivered by exactly one route in K (potentially different routes for the pickup and the delivery) and ( 2) routes in K can be scheduled at the cross-dock so that time constraints are respected. The overall algorithm is presented in 2. In Section 4.5.1, we present the dominance pre-processing we apply to routes (l. 1-4). The set partitioning (l. 6) is solved using a technique called branch-and-check, presented in Section 4.5.2. Its application to the VRPCD-RC is detailed in Section 4.5.3. The post-processing steps aimed at minimizing the volume transferred at the cross-dock in the obtained solution is presented in 4.5.3.

Dominance pre-processing

In the LNS, we stored both the routes and the legs in memory components. The goal of the later is to help improve the routes in K before solving the set partitioning. Algorithm 2: Perform SPS For that we identify non dominated legs. A pickup (resp. delivery) leg l i is said to be dominated by a leg l j if and only if: l i and l j serve the same set of requests, c j < c i and a j ≤ a i (resp. b j ≥ b i ) where c represents the cost of the leg, a represents the arrival time at the cross-dock ( resp. b represents the departure time from the cross-dock). Each route in K is improved by substituting its pickup and delivery leg by their non dominated equivalent.

Branch-and-check

We present branch-and-check [START_REF] Thorsteinsson | Branch-and-check: A hybrid framework integrating mixed integer programming and constraint logic programming[END_REF] using the following optimization problem:

M 1 : min c x (10) Ax ≤ b (11) H(x, y) (12)
x ∈ {0, 1} n (13)

y ∈ R m (14)
Assume that H(x, y) represents a set of constraints that have a limited impact on the LP relaxation and/or are difficult to efficiently model in a MIP, but that can be handled relatively easily by a CP solver. ( 10), ( 11) and ( 13) form a relaxation (M2) of (M1) that can be solved using branch-and-bound. The general principle of branch-and-check is the following. To solve (M1), a branch-and-bound is carried on (M2). Whenever an integral solution of (M2) is found, a CP solver is called to check constraints ( 12). If they are satisfied, the best solution found so far for (M1) is updated accordingly. Otherwise, this solution is rejected. In both cases the branch-and-bound process continues.

Application of branch-and-check to the VRPCD-RC

For the SPS in the VRPCD-RC, a classical set partitioning problem (SPP) is used as relaxation M2. For each request r ∈ R and each route k ∈ K, let λ + r,k (resp. λ - r,k ) be a binary constant which is equal to 1 if and only if this request is picked up (resp. delivered) by this route. And, for each route, let x r be a Boolean variable that indicate whether this route is selected. The SPP on routes is then :

min k∈K c k x k (15) k∈K λ + r,k x k = 1 ∀r ∈ R (16) k∈K λ - r,k x k = 1 ∀r ∈ R (17) x k ∈ {0, 1} ∀k ∈ K (18)
The objective ( 15) is to minimize the cost of the selected routes while constraints (16) (resp. ( 17)) ensure that each pickup point (resp. delivery point) is covered by exactly one routes.

A solution to the SPP on routes is a solution to the VRPCD-RC if and only if it is possible to find a schedule at the cross-dock such that time related constraints are respected. This is the problem presented in Section 4.4.1. In [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], on top of checking synchronization constraints, the set partitioning component was also assembling routes from legs (a process we called matching). Preliminary tests showed that, very often, the CP solver could not find any solution within reasonable run-times when asked to perform the matching with the extra resource constraints. Because the set partitioning idea proved very efficient in [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], we resort to keeping it but simplifying it to only do the scheduling part.

Improving legs pairing

The SPS approach is efficient to select the best routes, but it may still be possible to reduce the number of transfers at the cross-dock by improving the matching of legs to form routes in the best found solution. This can free processing capacity at the cross-dock and further help the search in subsequent steps. The solution obtained from the SPS is made of a set of pickup legs denoted Lp and a set of delivery legs denoted Ld . For each pickup leg l ∈ Lp , let T l be the set of delivery legs that deliver at least one request picked up by l. If a pickup leg l and a delivery leg l are matched together to create a route, there exist an associated unloading task o - ll , with a set of requests R - ll being unloaded, and a reloading task o + ll , with a set of requests R + ll being reloaded. These tasks have to be performed if and only if l and l are in the same route. Once again, the scheduling related concept used in CP can help with the modelling. Here, we each for each leg, exactly one of the previously described unloading (resp. reloading) task should be performed. This corresponds to alternative activities [START_REF] Beck | Scheduling alternative activities[END_REF] and can be modelled using alternative constraints (from IBM Corporation ( 2014)):

'An alternative constraint between an interval variable a and a set of interval variables b 1 , .

t l .Start ≥ t l .End ∀l ∈ Lp , l ∈ T l (25) o + ll .Start ≥ o - ll .End + δ r ∀l ∈ Lp , ∀l ∈ Ld s.t. R + ll = ∅ (26) o + ll .Start ≥ o - ll .End ∀l ∈ Lp , ∀l ∈ Ld s.t. R + ll = ∅ (27) o - ll .Start ≥ a l + δ u ∀l ∈ Lp , ∀l ∈ Ld s.t. R - ll = ∅ (28) o - ll .Start ≥ a l ∀l ∈ Lp , ∀l ∈ Ld s.t. R - ll = ∅ (29) o + ll .End ≤ b l ∀l ∈ Lp , ∀l ∈ Ld (30) l∈Lp isActive(t l , h) ≤ I ∀h ∈ [0, H] (31) l ∈L d isActive(t l , h) ≤ O ∀h ∈ [0, H] (32)
In this model, for each pickup leg l, t l is an interval variable that represents the associated unloading task that takes place at the cross dock. ( 19) minimizes the volume transferred at the cross-dock, while constraints (31). Alternative constraints (20) and ( 22) ensure that for each pickup leg l exactly one unloading task o ll is scheduled and that it is equal to t l . The same holds for delivery legs and reloading operations through variables t l and constraints ( 21) and (23). Constraints (24) ensure that the unloading operation associated with the matching of pickup leg l and the delivery leg l in the same vehicle is present if and only if the corresponding reloading operation is present as well. Constraints (25) ensure that all the reloading operations that depend on a pickup leg l start no earlier than the end of the unloading task associated with l. Constraints ( 26) and ( 27) ensure that when two legs are packed together, the delay between the two tasks respects the model presented in Section 3.2. Constraints ( 28) and ( 29) ensure that for each pickup leg, its corresponding unloading operation cannot start before the earliest feasible arrival time at the cross-dock. Constraints (30) ensure that for each delivery leg, its corresponding reloading operation is done by its latest feasible departure time. Constraints (32-31) account for the capacity constraints. A similar problem for the separated case can be formulated with an adaptation of (9). With the exception of the capacity constraints, the model is taken from [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF].

Computational experiments

The algorithm is coded in C++ and uses CPLEX and CP Optimizer from IBM ILOG Cplex Optimization Studio 12.6.1 as MIP solver and CP solver, respectively. The experiments were conducted under Linux using an Intel Xeon X7350 @ 2.93 GHz. Only one core is used both by our code and third party solvers. We consider instances proposed by [START_REF] Wen | Vehicle routing with crossdocking[END_REF], that range from 50 to 200 requests. They are based on real life data from a Danish logistics company. The termination criterion for all algorithms is based on time: 15 minutes for size 50 instances, 30 minutes for size 100 for instances, 60 minutes for size 150 instances and 120 minutes for size 200 instances. SPS time limit is 180 seconds. These time limits are between two and three times the run-times reported for the VRPCD in [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF]. As in [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], the number Φ of requests to remove in the repair phase of the LNS is drawn randomly in the interval [min(30, 10% of |R|), max(60, 20% of |R|)], acceptance criterion is descent.

Bound setting for the number of docks

To determine when limiting the number of docks start being a constraint in the VRPCD-RC, we post-processed the solutions obtained in [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD. To that end, for all the ten solutions found for each instance for the VRPCD, we solve optimization versions of the satisfaction scheduling problems introduced in Section 4.4.1. We take as objective: to minimize S in the shared case, and in the separated case, we only consider symmetric configurations where I = O and we minimize I. In Table 1, columns A correspond to the smallest value obtained after two hours of run-time for the CP solver, and columns B correspond to the worst value obtained after five minutes of post-process. For each instance, columns B corresponds to a threshold for the dock value, above which, the VRPCD-RC could be solved as a VRPCD (post-processing the solution for a limited amount of time to satisfy dock resource constraints). As such, in our experiments for the VRPCD-RC, we test dock values up to those reported in columns B. From this table, we can observe that dock resource constraints arise for dock values that corresponds to approximately 15% of the fleet size in the shared case and 10% of the fleet size in the separated case.

Parameters tuning

In this section we evaluate and adjust several parameters. We first give some insights on the impact of parameters. We start with the time limit for the CP solver in LNS feasibility tests (called CP time limit), then we report the success rate of heuristics. After, we present the influence of the SPS frequency on the quality of solutions. Eventually, we compare the performance of four possible configurations: with CP solver tests/with heuristic tests in LNS, with/without SPS. For tuning, we use instances 50b, 100b, 150b, 200b, and we consider the shared cross-dock configuration case. Parameter tuning has been performed taking an overall time limit for LNS or LNS+SPS (called LNS time limit) and comparing the final average cost value over these representative instances. Table 1: Dock value obtained when post-processing for each instance all of ten solutions of [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF]. Columns A refer to the best solution (min dock use) obtained after two hours, while columns B correspond to the worst value (max dock use) obtained after five minutes of post-processing.

CP solver time limit in feasibility tests

When checking the feasibility with respect to dock resource constraints, we need to set a time limit after which the CP solver will stop searching and declare the insertion unfeasible. This so-called CP time limit prevents spending a large amount of time checking the feasibility of a single insertion. Four values have been tested and the final result was that a CP-time limit of 0.01 seconds yields better results. To understand this observation, Table 2 provides some insights on the average percentage of unsolved feasibility problems for each evaluated CP time limit.

CP Time Limit (s) 1 0.1 0.01 0.001 Run-time 1 0.33 0.12 0.10 Time limit hit (%) 0.5 12.9 20.1 27.3

Table 2: Comparison of four different time limits for the CP solver, when used as feasibility test in the shared case. Run-time is normalized with 1 representing the run-time for a CP solver time limit of 1s. Time limit hit represents the percentage of calls for which the CP solver could not find an answer within the time limit. Figures reported for one thousand LNS iterations, two runs were performed in each case.

Insights on heuristics performance

As mentioned in Section 4.4.2, using CP versus using heuristics for feasibility tests is a quality/run-time trade-off. For each instance in the training set, we performed two runs with LNS, with the stop-criterion set to one thousand iterations. We count how many time insertions were reported feasible by heuristics and how many times they were reported feasible by the CP solver with a time limit of 1 second (according to Table 2, with this CP time limit we can consider that the CP solver give an accurate answer most of the time). On average, 71.1% of feasible insertions are reported feasible by heuristics. Performing one thousand LNS iterations with feasibility tests performed by heuristics takes only 18.2% of the time taken with CP solver tests (with 0.01s as CP time limit).

Insights on SPS frequency

In [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], the SPM was solved twenty times per run (every thousand iterations with a stop-criterion of twenty thousands iterations). In the VRPCD-RC, because feasibility tests are computationally intensive, the stop-criterion is based on time. The number of iterations performed within the LNS time limit depends not only on the size of the instance but also on the dock value. As such, we propose a SPS-criterion based on time. In Table 3 andTable 4

Performance comparison

On Fig. 3 we present the convergence curves of LNS and LNS+SPS over time for both CP solver and heuristics feasibility tests. From this graph, we can observe two things. First, as in the VRPCD, periodically solving a set partitioning based problem significantly improves performance compared to LNS alone: -6.96% for heuristic test and -7.92% for CP solver test on average. Second, the best performing method is LNS+SPS with heuristic feasibility tests, as it finds solutions that are 1.19% better on average than LNS+SPS with CP solver test. In the rest, we thus use LNS+SPS with heuristic feasibility test.

Results

Table 5 presents a synthesis of the results in both the shared case and separated cases respectively (detailed results are presented in Appendix A). A dock value of 0 corresponds to the case where no transfer is allowed, i.e. vehicles have to go back to the cross-dock between pickup and delivery legs but cannot transfer any items. We present this value mostly for comparison purposes. The algorithm shows relatively good performance in terms of stability: the difference between the average value and the best value of the five runs for each instance and each dock value is 0.6% on average and at most 1.4% in the shared case and 0.6% on average and at most 3.2% in the separated case.

Regarding routing costs: integrating dock resource constraints implies an increase in cost. Comparing the two systems, the shared case costs are slightly smaller than separated costs (e.g, for instance 150b, 5.85% for 2 shared docks compared to 6.05 as shown on Fig. 4). These differences increases with the number of docks.

Conclusion

This paper presents an adaptation of the matheuristic proposed in [START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF], which is based on large neighborhood search and periodic calls to a set partitioning based problem, to solve an extension of the VRPCD that includes resource synchronization constraints at the cross-dock. To deal with these constraints, scheduling heuristics and a CP model have been used as feasibility tests for insertions in LNS. In this case, experiments have showed that heuristics are the most efficient compromise in terms of run-time versus solution-quality. The proposed method has 6 and7 been tested on instances from the literature, where it shows an increase in routing costs with the decrease in cross-dock capacity.

Adding resource synchronization constraints has made the problem computationally very challenging: several decision had to be taken (heuristics as feasibility tests, simpler set partitioning component) to limit the run-time. This shows that integrated problems are challenging even for heuristics. From the application point of view, several applications and extension of the proposed method could be envisioned in the field of electric vehicle routing [START_REF] Schneider | The electric vehicle-routing problem with time windows and recharging stations[END_REF][START_REF] Enderer | Integrating dock-door assignment and vehicle routing with cross-docking[END_REF] and city logistics [START_REF] Grangier | An adaptive large neighborhood search for the two-echelon multiple-trip vehicle routing problem with satellite synchronization[END_REF].

Fig. 1 :

 1 Fig. 1: Time chart for vehicle unloading and reloading at the cross-dock in the VRPCD-RC

Result:

  The best found solution s 1 Pool of legs L := ∅ 2 Pool of routes K := ∅ 3 Generate an initial solution s 4 s := s 5 while stop-criterion not met do 6 s := s 7 Destroy quantity: select a number Φ of requests to remove from s 8 Operator selection: select a destruction operator M -and a repair operator M + 9 Destruction : apply M -to remove Φ requests from s , and put them in the requests bank of s 10 Repair: apply M + to reinsert the requests in the requests bank in s 11 if acceptance criteria is met then 12 s := s 13 end 14 if cost of s is better than cost of s then 15 s := s 16 end 17 Add legs of s to L 18 Add routes of s to R 19 if set partitioning and scheduling condition is met then 20 Perform SPS(L, K, s , s );

Fig. 2 :

 2 Fig. 2: Logical flow-chart of feasibility tests for the VRPCD-RC

Fig. 3 :

 3 Fig. 3: Comparison of the evolution of the average solution quality over time (in percentage of LNS time limit) for four different configurations: with CP solver/with heuristics tests, with/without SPS. The results aggregate 5 runs for each dock value in the shared case for instances in the training set. They have been normalized, first by instance then by method, with 100 representing the cost at the end of LNS with CP solver tests (LNS-CP)

  Fig.4: Influence of the dock value for instance 150b. Five runs were performed for each dock value. The y-axis represents the average gap with respect to the average value reported in[START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF]

  Input: pool of legs L, pool of routes K, solution s , s from Algorithm 1 1 Remove dominated legs from L 2 for each route k ∈ K do

	3	Replace, if possible, its pickup leg and/or its delivery leg with a non-dominated
		equivalent in L
	4 end
	5 Remove dominated routes from K
	6 Solve SPS with all routes in K
	7 if a new best solution has been found then
	8	Improve the matching of legs (as in Section 4.5)
	9	Update s
	10 end
	11 if Set partitioning was not solved to proven optimality then
	12	Clear K
	13 end

  . . , b n models an exclusive alternative between b 1 , . . . , b n . If interval a is present, then exactly one of intervals b 1 , . . . , b n is present and a starts and ends together with this specific interval. Interval a is absent if and only if all intervals in b 1 , . . . , b n are absent.'

	In the shared case, the post-processing is performed by solving the following CP
	model:
	min	o -ll .IsP resent ×	q r	(19)
	(l,l )∈ Lp× Ld	r∈R -l,l
	Alternative(t l , {o -ll ; ∀l ∈ Ld })	∀l ∈ Lp (20)
	Alternative(t l , {o + ll ; ∀l ∈ Lp })	∀l ∈ Ld (21)
		o -ll .IsOptional ← T rue	∀l ∈ Lp , ∀l ∈ Ld (22)
		o + ll .IsOptional ← T rue	∀l ∈ Lp , ∀l ∈ Ld (23)
	o -ll .IsP resent ⇐⇒ o + ll .IsP resent	∀l ∈ Lp , ∀l ∈ Ld (24)

  , we compare three different settings for the SPS frequency : 10, 20 and 40 calls within the LNS time limit. Accordingly, the number of calls per run is set to 20 for heuristic test and 10 for CP solver test.

	Number of calls	10	20	40
	Average gap (%)	0	-0.70 -0.29

Table 3 :

 3 Comparison of the impact of the number of SPS calls per run on the quality of the solution for heuristic feasibility tests. Five runs were performed for each dock value for each instance in the training set. 10 calls is taken as reference for the gap

	Number of calls	10	20	40
	Average gap (%)	0	+0.40 +0.68

Table 4 :

 4 Comparison of the impact of the number of SPS calls per run on the quality of the solution for CP solver tests. Five runs were performed for each dock value for each instance in the training set. 10 calls is taken as reference for the gap

Table 5 :

 5 Gaps (in %) in the shared and separated cross-dock configuration cases with respect to the average and best values reported in[START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD. A * indicates that not all instances had to be tested on such cross-dock values, as the best non dock-resource constrained solution did not require that many cross-docks. Detailed results per instances are presented in Tables

Table 7 :

 7 Average values and best solution found in the separated cross-dock configuration case; LNS+SPS was run five times for each instance. Columns Gap refer to the gap to average values and best solutions reported in[START_REF] Grangier | A matheuristic based on large neighborhood search for the vehicle routing problem with cross-docking[END_REF] for the VRPCD

	Instance Dock	Average Value Gap (%)	Value	Best Gap (%)	Instance Dock	Average Value Gap (%)	Value	Best Gap (%)
	50a	0 1	6882.4 6682.72	6.5 3.4	6871.93 6628.05	6.4 2.7	0 1	21036.12 20827.96	7.9 6.8	21019.3 20703.5	8.1 6.4
		2	6530.81	1.0	6471.48	0.2	2
	50b	1 0	7609.07 8027.87	2.4 8.1	7541.33 8021.81	3.0 9.6	150e
		2	7481.99	0.7	7470.78	2.0
	50c	0 1	7857.77 7505.34	7.3 2.5	7827.67 7473.94	7.1 2.2
		2	7449.88	1.8	7397.38	1.2
		3	7338.88	0.3	7330.08	0.3
	50d	0 1	7760.32 7272.57	10.2 3.3	7760.11 7206.48	10.4 2.5
		2	7108.91	1.0	7076.13	0.7
		0	8157.77	9.1	8156.94	9.4
	50e	1	7843.08	4.9	7772.39	4.3
		2	7616.54	1.8	7554.75	1.4
		3	7537.85	0.8	7473.56	0.3
		0	15636.26	8.8	15628.7	8.9
	100b	1 2	15322.86 15181.04	6.6 5.6	15234.2 15073.8	6.2 5.0
		3	14929.18	3.8	14810.3	3.2
		4	14700.96	2.3	14620.3	1.9
		0	14915.92	7.9	14915.6	8.2
	100c	1 2	14654.16 14456.6	6.0 4.5	14611.5 14395.0	6.0 4.4
		3	14353.9	3.8	14188.2	2.9
		4	14145.78	2.3	14081.3	2.2
		0	14860.46	9.3	14832.5	9.2
	100d	1 2	14424.74 14316.7	6.1 5.3	14338.9 14207.4	5.6 4.6
		3	14070.18	3.5	13976.2	2.9
		4	13866.24	2.0	13774.7	1.5
		0	15095.92	8.1	15091.5	8.2
		1	14819.26	6.2	14733.8	5.7
	100e	2	14783.82	5.9	14622.8	4.9
		3	14515.0	4.0	14384.7	3.2
		4	14357.26	2.9	14258.4	2.3
		5	14281.0	2.3	14085.3	1.0
		0	20859.22	7.5	20807.7	7.5
		1	20534.8	5.8	20406.4	5.4
		2	20505.68	5.7	20342.5	5.1
	150a	3	20248.54	4.4	20090.1	3.8
		4	20110.14	3.7	19988.6	3.3
		5	19822.08	2.2	19740.0	2.0
		6	19725.04	1.7	19629.7	1.4
		7	19603.78	1.0	19541.6	0.9
		0	22272.64	7.7	22236.5	8.0
		1	22018.68	6.5	21934.5	6.6
		2	21882.36	5.9	21825.0	6.0
	150b	3	21651.48	4.7	21585.3	4.9
		4	21432.16	3.7	21360.2	3.8
		5	21227.66	2.7	21097.6	2.5
		6	21065.9	1.9	20986.2	2.0
		7	20963.44	1.4	20899.9	1.5
		0	21313.88	7.8	21295.7	8.0
		1	21095.74	6.7	20959.6	6.2
		2	20945.68	5.9	20854.4	5.7
	150c	3	20869.2	5.5	20642.1	4.6
		4	20685.7	4.6	20559.7	4.2
		5	20475.68	3.6	20369.1	3.3
		6	20208.76	2.2	20155.8	2.2
		7	20153.68	1.9	20100.9	1.9
		0	21962.58	7.9	21951.3	8.0
		1	21735.46	6.8	21536.0	6.0
		2	21609.34	6.2	21442.6	5.5
	150d	3	21448.08	5.4	21368.0	5.2
		4	21504.76	5.6	21289.8	4.8
		5	21040.52	3.4	20874.7	2.7
		6	20806.38	2.2	20642.2	1.6
		7	20715.78	1.8	20608.5	1.4

Acknowledgments

This work was partially supported by the Canadian Natural Science and Engineering Research Council (RGPIN-2015-04696 ) and by the Fonds de recherche du Québec -Nature et technologies through its Team research Program.

Appendix A Detailed results

Table 6 and Table 7 present the detailed results obtained for each instance. These results highlight the increased complexity induced by integrating dock resource constraints, as there exists for each instance, at maximum dock value, a solution with 0% gap (see Section 5.1). Nevertheless LNS+SPS finds solutions that are 1.6% more expensive in the shared case and 1.5% more expensive in the separated case, which remains satisfactory.


Gao-Ren Li
email: ligaoren@mail.sysu.edu.cn

Dun-Lin Qu

Zi-Long Wang

Cheng-Yong Su

Ye-Xiang Tong

Laurent Arurault

Ceria-terbia solid solution nanobelts with high catalytic activities for CO oxidationw

Keywords: Catalytic activity measurements, characterization, EDS pattern, XPS spectrum, SEM images, TG curves

Ceria-terbia solid solution nanobelts were prepared by an electrochemical route and tested as catalysts of high activity for CO oxidation.

Air pollution has become one of the major global problems. Three-way catalysts (TWC) for automotive exhaust treatment, which can transform three environmentally hazardous gases, namely, CO, NOx, and hydrocarbons to environmentally benign or less hazardous form, is therefore a topic of intensive research. [1][2][3][4][5] Ceria-based materials are very important components of promoters for automotive TWCs. [6][7][8][9] It is well known that rare earth ion-doped ceria can increase the ionic conductivities, 10 and therefore may generally enhance the catalytic activity of ceria by making oxygen from the bulk more accessible for surface reactions. 11 As is well known, Ce 4+ and Ce 3+ always coexist in CeO 2 . 12 The ionic radius of Ce 3+ ion (1.14 A ˚) is larger than that of Ce 4+ ion (0.97 A ˚) and so the surrounding fluorite lattice will be distorted when Ce 3+ vacancies are formed in CeO 2 . When different +4 cations are doped into the CeO 2 lattice, the oxygen diffusivity within the lattices will be influenced. The energetic cost of these distortions in CeO 2 can be decreased by formation of solid solutions with 4+ ions of smaller size than Ce 4+ ions, such as Zr 4+ , Hf 4+ , Pr 4+ , Tb 4+ . [11][12][13][14][15] In addition, the ion mobility inside the CeO 2 lattice will be changed, resulting in the formation of a defective fluorite-structured solid solution. Such modifications in the defect structure of CeO 2 will confer new properties to the catalyst such as high catalytic activity. 16 The use of variable valence dopants in CeO 2 has attracted much attention. [17][18] Herein, special emphasis is given on the influence of Tb 4+ doping on the catalytic performance of ceria-based mixed oxides. The element Tb can exhibit variable valence states (3+ and 4+), which may have some special impact on the final product.

Constructing novel nanostructures to obtain enhanced properties appears a viable route in designing and selecting catalysts and catalyst promoters for specific catalytic reactions. 19,20 CeO 2 -based composites with nanostructures usually have high surface areas, improved sintering properties and high oxygen storage capacities, which may highly enhance the catalytic performance. To the best of our knowledge, there have been no reports on the successful synthesis of onedimensional ceria-based mixed oxide nanobelts. Herein, detached Ce 1Àx Tb x O 2Àd nanobelts are first synthesized by electrochemical deposition at room temperature and usual atmospheric pressure without any surface-capping agent. The Ce 1Àx Tb x O 2Àd nanobelts as a new family of CeO 2 -based oxide nanostructures are desirable because of their unique physical and chemical properties. 11 The electrochemical deposition is a simple and low-cost route to this material.

In this experiment a simple three-electrode cell was used. A highly pure Pt foil (99.99 wt%, 0.25 cm 2 ) was used as the auxiliary electrode. A saturated calomel electrode (SCE) was used as the reference electrode that was connected to the cell with a double salt bridge system. The electrochemical deposition was carried out in a solution of 0.01 M Ce(NO 3 ) 3 + 0.001 M Tb(NO 3 ) 3 + 0.1 M NH 4 NO 3 with a current density of 4.0 mA cm À2 in a configured glass cell at room temperature, in which a Cu plate (99.99 wt%, 0.25 cm 2 ) served as the substrate. The oxygen was dissolved in the deposition solution through bubbling in the electrolyte. Before electrodeposition, the Cu substrate was cleaned ultrasonically in 0.1 M HCl, distilled water, and acetone and then rinsed in distilled water again. The products were characterized by X-ray diffractometry (D/MAX 2200 VPC with Cu-Ka radiation), and an Oxford Instruments INCA energy-dispersive spectrometer (EDS) was employed to analyze chemical composition. Microstructures of the deposits were characterized by field emission scanning electron microscopy (FE-SEM; JSM-6330-F) and X-ray photoelectron spectroscopy (XPS, ESCALAB 250) was used to assess the chemical state and surface composition of the deposits.

SEM images of Ce 1Àx Tb x O 2Àd composites with different magnifications are shown in Fig. 1(a) and (b), which reveal that the deposit consists of a large quantity of nanobelts with length of about 0.5-2.0 mm, width of 100-250 nm and thickness of only about 30 nm. The microstructural details of Ce 1Àx Tb x O 2Àd nanobelts were further investigated by TEM as shown in Fig. 1(c) and (d). High-resolution TEM (HRTEM) images of individual nanobelts clearly display lattice fringes for nanocrystals, indicating that these Ce 1Àx Tb x O 2Àd nanobelts possess high crystallinity. The measured interplanar spacings from HRTEM are about 0.31 and 0.28 nm, which are in good agreement with (111) and (200) facet distances of the CeO 2 phase, respectively. The SAED pattern (inset in Fig. 1(d)) indicates that the prepared nanobelts have a polycrystalline structure. It is well known that a small amount of Ce 3+ coexists at the surface of CeO 2 . 21 Two sets of spin-orbit multiplets corresponding to 3d 3/2 and 3d 5/2 are present as peaks labelled u and v, respectively in Fig. 2(a). The peaks centered at 899.0 eV (u 1 ), 908.3 eV (u 2 ) and 917.1 eV (u 3 ) eV can be attributed to the Ce 4+ contribution while the peaks centered at 884.2 eV (v 1 ) and 902.2 eV (v 2 ) can be attributed to the Ce 3+ contribution. Thus, the surface of the present samples contains Ce both in the 4+ and the 3+ states. Fig. 2(b) shows the Tb 4d spectrum of a Ce 1Àx Tb x O 2Àd nanobelt sample. As the kinetic energies of the photoemitted electrons are very low, Tb 3d core level analysis is very difficult by a conventional laboratory spectrometer (with either Mg or Al KR radiation). Thus the most intense core level (Tb 4d) is used for analysis. There is an agreement that Tb 3+ gives a signal below 150 eV and Tb 4+ is related to features well above 150 eV. [22][23][24] The spectrum in Fig. 2(b) obviously shows a pronounced shoulder below 150 eV, a maximum at 152.3 eV and a tailing toward 160 eV, indicating the presence of more than one oxidation state, i.e., the 3+ and the 4+ oxidation states. Therefore, the XPS results further demonstrated Ce 1Àx Tb x O 2Àd composites were successfully synthesized.

Fig. 3(a) shows the XRD pattern of as-prepared Ce 1Àx Tb x O 2Àd (x = 6.5 at.%) nanobelts. Peaks corresponding to CeO 2 (111), ( 200), ( 220) and (311) planes were observed, indicating a face-centered cubic-phase CeO 2 (JCPDS 34-0394). The diffraction peaks of TbO 2 or Tb 2 O 3 were not observed, indicating Tb has entered into the CeO 2 lattice. The calculated cell parameter (a) is equal to 0.5378 nm, a little smaller than that of bulk CeO 2 (0.5411 nm). This may be due to a lattice constriction effect resulting from Tb 4+ ions as Tb 4+ has a smaller ionic radius than that of Ce 4+ (the ionic radii values for Ce 4+ , Tb 4+ , Ce 3+ and Tb 3+ are 0.97, 0.88, 1.14 and 1.04 A ˚, respectively). 25 The lattice constriction is additional evidence for the formation of solid solutions. The broadening of the reflections of Ce 1Àx Tb x O 2Àd nanobelts in Fig. 3(a) can be attributed to the smaller nanocrystallites. The average crystallite sizes of Ce 1Àx Tb x O 2Àd nanobelts were calculated from X-ray line broadening of the (111) reflection using Scherrer's equation (D = Kl/(bcos y), where l is the wavelength of the X-ray radiation, K is a constant, y is the diffraction angle, and b is the full width at half-maximum). The average crystallite size of Ce 1Àx Tb x O 2Àd nanobelts was determined to be about 5.5 nm. The peak shift of Ce x Tb 1Àx O 2Àd (x = 6.5 at.%) nanobelts compared with pure CeO 2 can be estimated as about 0.41. The XRD pattern of Ce 1Àx Tb x O 2Àd nanobelts calcined at 1100 K are shown in Fig. 3(b). The same monophasic composition is retained, and no phase segregation is observed.

The catalytic activity of the Ce 1Àx Tb x O 2Àd (x = 6.5 at.%) nanobelts was investigated by performing the CO conversion reaction. Fig. 4 shows the catalytic activities of Ce 1Àx Tb x O 2Àd nanobelts, nanoparticles and nanosheets (all x = 6.5 at.%) as a function of reaction temperature (the surface areas of Ce x Tb 1Àx O 2Àd nanobelts, nanosheets and nanoparticles are about 115, 72 and 23 m 2 g À1 , respectively). It can be clearly observed the CO conversion increases with increasing reaction temperature for all samples. For Ce 1Àx Tb x O 2Àd nanobelts as catalyst, a 96% CO conversion is achieved at about 913 K. However, only 43 and 57% CO conversions are obtained at Since the components of Ce 1Àx Tb x O 2Àd deposits were kept the same, the above different catalytic activities must originate from the different structures that have different oxygen storage capacities (OSCs). As is well known, the value of OSC is a crucial parameter for catalytic activity. Herein, the OSC properties of Ce 1Àx Tb x O 2Àd nanobelts, nanoparticles and nanosheets were tested by oxygen release characteristics of the samples calcined at different temperatures under dry air atmosphere in the temperature range 450-1100 K. The change of weight of the sample was monitored by thermogravimetry (TG) under cyclic heat treatments in flowing air. Ce 1Àx Tb x O 2Àd samples were subjected to consecutive cycles of heating and cooling, and typical TG curves of various samples were shown in Fig. S5 (ESIw). In the above results, a large decrease of weight was observed in the first heat treatment (A -B) and this can be attributed to the release of both water molecules and oxygen from Ce 1Àx Tb x O 2Àd solid solutions. In the cooling back stage (B -C), a partial recovery of the weight of Ce 1Àx Tb x O 2Àd sample was observed. In the second heating cycle (C -B), a small decrease of weight was also observed corresponding to the potential oxygen release capacity of the powders in ordinary air atmosphere. 12 The Ce 1Àx Tb x O 2Àd nanobelts show the highest weight loss, which is followed by Ce 1Àx Tb x O 2Àd nanosheets and nanoparticles. Accordingly the Ce 1Àx Tb x O 2Àd nanobelts have the highest OSC among various Ce 1Àx Tb x O 2Àd nanostructures. These results are consistent with the catalytic activity order of CO oxidation, where Ce 1Àx Tb x O 2Àd nanobelts showed better activity than that of Ce 1Àx Tb x O 2Àd nanosheets, followed by Ce 1Àx Tb x O 2Àd nanoparticles.

In summary, a novel design of ceria-terbia solid solution nanobelts led to catalysts of high activity for CO oxidation.

Notes and references

  A representative EDS pattern of Ce 1Àx Tb x O 2Àd nanobelts is shown in Fig. S1 (ESIw) which revealed O, Ce and Tb to be present. The content of Tb in Ce 1Àx Tb x O 2Àd composites was determined to be about 6.5 at.%, and the Ce content is about 93.5 at.%. XPS analyses of Ce 1Àx Tb x O 2Àd nanobelts were carried out and the O 1s core level XPS profile is shown in Fig. S2 (ESIw). The peak centered at 530.1 eV corresponds to the O 2À contribution, and can be attributed to the lattice oxygen associated with Ce 1Àx Tb x O 2Àd . Another peak at a slightly lower binding energy side can be attributed to absorbed oxygen species (presumably from absorbed water).

Fig. 1

 1 Fig. 1 SEM images of (a) low magnification and (b) high magnification; (c) TEM image; (d) HRTEM image; and SAED (inset) of Ce 1Àx Tb x O 2Àd nanobelts prepared in solution of 0.01M Ce(NO 3 ) 3 + 0.001 M Tb(NO 3 ) 3 + 0.1M NH 4 NO 3 with current density of 4.0 mA cm À2 .

Fig. 3

 3 Fig. 3 XRD patterns of the obtained samples: (a) as-prepared Ce 1Àx Tb x O 2Àd nanobelts; (b) Ce 1Àx T Àx O 2Àd nanobelts calcined at 1100 K.

  Ce 1Àx Tb x O 2Àd nanobelts (x = 6.5 at.%) have been synthesized by electrochemical deposition at room temperature and ordinary atmospheric pressure without any surface capping agent. The electrochemical deposition route shows a facile and low-cost route. The results of XRD and TEM show Ce 1Àx Tb x O 2Àd nanobelts are well crystallized and have a cubic structure. A strong shape/crystal structure effect of Ce 1Àx Tb x O 2Àd composites for the catalytic activity has been identified. The Ce 1Àx Tb x O 2Àd nanobelts are most active for the CO conversion reaction compared with Ce 1Àx Tb x O 2Àd nanoparticles and nanosheets. The high catalytic activity of Ce 1Àx Tb x O 2Àd nanobelts can be attributed to the special nanostructures with high OSC. These findings can be extended to new studies of other ceria-based mixed oxides with high surface areas. This work was supported by NSFC (20603048, 20873184 and 90923008) and Guangdong Province (2008B010600040 and 9251027501000002).

Fig. 4

 4 Fig. 4 Conversion of CO over (a) Ce 1Àx Tb x O 2Àd nanobelts, (b) Ce 1Àx Tb x O 2Àd nanosheets and (c) Ce 1Àx Tb x O 2Àd nanoparticles as a function of reaction temperature.