Pawan Kumar

Harshal P Mungse

Stéphane Cordier

Rabah Boukherroub

Om P Khatri

Suman L Jain
email: sumanjain@iip.res.in

Hexamolybdenum Clusters Supported on Graphene Oxide: Visible-Light Induced Photocatalytic Reduction of Carbon Dioxide into Methanol

Keywords: Graphene oxide, Mo cluster compound, photocatalysis, carbon dioxide reduction, methanol

Hexamolybdenum (Mo 6 ) cluster-based compounds namely Cs 2 Mo 6 Br i 8 Br a 6 and (TBA) 2 Mo 6 Br i 8 Br a 6 (TBA = tetrabutylammonium ) were immobilized on graphene oxide (GO) nanosheets by taking advantage of the high lability of the apical bromide ions with oxygen-functionalities of GO nanosheets. The loading of Mo 6 clusters on GO nanosheets was probed by Fourier-transform infrared (FTIR) spectroscopy, X-ray photoelectron spectroscopy (XPS), high resolution transmission electron microscopy (HRTEM) and elemental mapping analyses. The developed GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x composites were then used as heterogeneous photocatalysts for the reduction of CO 2 under visible light irradiation. After 24 h visible light illumination, the yield of methanol was found to be 1644 and 1294 µmol.g -1 cat for GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x , respectively. The quantum yields of methanol by using GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x as catalysts with reference to Mo 6 cluster units presented in 0.1g amount of catalyst were found to be 0.015 and 0.011, respectively. The role of immobilized Mo 6 clusters-based compounds on GO nanosheets is discussed to understand the photocatalytic mechanism of CO 2 reduction into methanol.

Introduction

Graphene oxide (GO), an oxidized form of graphene, conventionally served as a most established precursor for the preparation of reduced graphene oxide, chemically functionalized graphene and graphene-based composites [1][2]. Owing to its rich surface chemistry, high specific surface area and various structural defects, GO has attracted large attention not only as a catalytic material, but also as an excellent nanostructured support for heterogenization of various catalysts for efficient chemical transformations [3][4]. The first report on the catalytic activity of GO was revealed by Bielawski et al. in 2010 for the oxidation of various alcohols and the hydration of various alkynes [5]. Then after, the catalytic performance of GO was expanded for various chemical transformations, including oxidation of thiols, sulphides and olefins, aza-Michael addition of amines to activated alkenes, acetalization of aldehydes, ring-opening polymerization of various cyclic lactones and lactams, Friedel-Crafts addition of indoles to α,β-unsaturated ketones, oxidative coupling of amines to imines, oxidative amidation of aromatic aldehydes, etc [6][7][8][9][10][11]. The catalytic properties of GO have mainly originated from the presence of various oxygen functionalities and unpaired electrons in the GO scaffold and its very high surface area [2][3]. Recently, Riedo et al have demonstrated that GO is a metastable material, whose structure and chemical functionalities evolve at room temperature [12]. Furthermore, an increase of temperature induces gradual reduction of oxygen functionalities in the GO [13]. The ease of GO aging and reduction during catalytic reactions under mild to high temperature raised the issue of recyclability of GO for chemical transformations [14].

Moreover, the high specific surface area, excellent dispersibility in water/ polar solvents, presence of plenty of oxygen containing functional groups and several structural defects in GO scaffold provide ample opportunities to serve as a nano-structured support for heterogenization of catalytically active metal complexes, organic functionalities, and transition metal clusters for efficient chemical transformations [15][16][17][18][19].

In recent years, several methodologies have been established for immobilization of various nano-materials on the GO nanosheets, which include the use of bi-functional linkers, light-induced immobilization, co-reduction of metal salts in the presence of GO, etc [2,[20][21][22]. In most of these approaches, the oxygen functionalities and structural defects in the GO scaffold function as anchoring sites for the immobilization of various types of nanomaterials.

Molybdenum octahedral (Mo 6 ) clusters with the general formula A 2 Mo 6 X i 8 X a 6 (A = Cs + , n-(C 4 H 9 ) 4 N + and X = Br -) exhibit photochemical and redox properties due to delocalization of valence electrons on the metal centers [23][24][25]. The Mo 6 clusters are face-capped by eight inner ligands (Br i ) through Mo-Br i bonds having covalent dominant character, and are additionally bonded to six apical ligands (Br a ) through Mo-Br a linkage having strong ionic character. The dissolution of these clusters affords [Mo 6 Br i 8 Br a 6 ] 2-anionic units, where the Mo 6 Br i 8 core constitutes a rigid block that can be functionalized by exchange / replacement of Br a located at the apical position of the cluster. Owing to their remarkable structural and photo-physical properties, A 2 Mo 6 X i 8 X a 6 cluster-based compounds provide immense potential for nanoarchitectonic and photocatalytic applications [26]. Recently, Boukherroub et al. have demonstrated that nanocomposites of Mo 6 clusters, gold nanoparticles and GO exhibited very high photocatalytic activity for the degradation of rhodamine B under visible light irradiation. In these photocatalysts, the Mo 6 cluster absorbs the incident photons and initiates the photocatalytic process [24]. The Mo 6 clusters as homogeneous catalyst could also afford photo-reduction of CO 2 into methanol induced by visible light. However, the need of sacrificial donor and non-recyclability of these clusters make the catalytic process expensive [27].

Due to the rising concentration of CO 2 and the depletion of fossil fuel reserves, there is considerable interest in the development of alternative energy sources [28][START_REF] Aresta | Carbon dioxide reduction and use as a chemical feedstock[END_REF][START_REF] Indrakanti | [END_REF]. Conversion of CO 2 into hydrocarbon fuels utilizing the abundant solar energy can provide a viable solution for both problems simultaneously [31][START_REF] Halmann | Greenhouse Gas Carbon Dioxide Mitigation Science and Technology[END_REF]. In this regard, a number of heterogeneous catalysts including metal oxides, sulphides, nitrides and oxynitrides have been investigated for CO 2 reduction, albeit lower quantum yields and poor selectivities are the major drawbacks [START_REF] Izumi | [END_REF][34][35]. Molecular catalysis, particularly, based on transition metal complexes, such as rhenium(I) bipyridine, ruthenium(II) polypyridine carbonyl, cobalt (II) trisbipyridine, and cobalt (III) macrocycles in combination with a photosensitizer have been considered to be high performance photocatalysts for CO 2 reduction with relatively high quantum yield and high selectivity of products [36][37][38][39]. However, the homogeneous nature and non-recycling ability of these catalysts make them impractical from economical as well as environmental viewpoints. This problem can be overcome by anchoring these complexes to photoactive supports, which not only provides facile recovery of the catalyst, but also photoactive support and complex may work synergistically for better electron transfer to CO 2 [40][41]. In a recent report, we have demonstrated the photocatalytic activity of Mo 6 clusters for the photoreduction of CO 2 using triethylamine as the sacrificial agent [27]. However, non-recycling ability of the photocatalyst and need of sacrificial donor (triethylamine) are the major drawbacks of the developed method. Herein, Cs into methanol (Figure 1). examined on high resolution transmission electron microscope (HRTEM) by drop casting their ethanol solutions on the TEM grids. Field emission scanning electron microscopy (FESEM) analyses were carried out using an FEI Quanta 200 F. The elemental distribution in the developed composite materials was recorded by using energy dispersive spectrometer (EDS) coupled with FESEM. UV/Vis spectra were measured on a Hitachi U-2900 UV/Vis spectrometer with a path length of 1 cm using quartz cuvettes.

Experimental

Photocatalytic Reduction of CO 2

Photocatalytic reduction of CO 2 was carried out using GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x composites as the photocatalysts. The experiments were performed in borosil cylindrical vessel (100 mL) of 5 cm diameter. Photo-irradiation was carried out under visible light by using 20 W white cold LED flood light (model no. HP-FL-20W-F-Hope LED Opto-Electric Co., Ltd.). The intensity of the light at vessel was measured by intensity meter and was found to be 75 Wm -2 . The vessel was charged initially with water (10 mL) and DMF (40 mL), and then the solution was degassed by continuous purging of nitrogen for 15 min. Subsequently, CO 2 was bubbled through the solution for at least 30 min to saturate the solution, then GO-Cs 2 Mo 6 Br i 8 Br a x or GO-(TBA) 2 Mo 6 Br i 8 Br a x catalyst (100 mg) was added to the above-described solution. The vessel was tightly closed during the photocatalytic reaction and stirred continuously by a magnetic stirring bar to prevent sedimentation of the catalyst. Samples were collected every 2 h by using a needle, and the catalyst was removed by filtration (0.2 µm PTFE, 13 mm diameter). Quantitative determination of reaction product was achieved by using GC-FID (Varian CP3800 by using CP Sil 24CB LOW BLEED/MS column, flow rate: 0.5 mL.min -1 , injector temperature: 250 °C, FID detector temperature: 275 °C). A calibration curve was generated for quantification and for confirmation of linear response of the GC-FID system. Further confirmation of methanol production was done by HPLC (Shimadzu UFLC, by using Oyster BDS Premium C18 250X4.6mm, 5µm column, mobile phase acetonitrile: acetone (65:35), flow rate: 0.5 mL/min at 205 nm).

Blank reactions were conducted to ensure that methanol production was due to the photo-reduction of CO 2 and to eliminate surrounding interference. For the determination of gaseous products, 20 µL sample was withdrawn and analyzed with the help of GC-TCD and GC-FID (Agilent 7890A GC system) equipped with RGA (refinery gas analyzer) capillary column; flow rate (H 2 : 35 mL/min, air: 350 mL/min, makeup flow: 27 mL/min, for TCD reference flow: 45 mL/min, Helium flow: 2mL/min), injector temperature: 220°C, TCD detector temperature and FID detector temperature: 220°C.

Results and discussion

Synthesis and characterization of catalyst

The chemical and structural features of GO and its composites i. the quantitative estimation along with chemical changes in the GO and Mo 6 clusters (Figures 3 and4). The high-resolution C 1s spectrum of GO exhibited an overlapped double peak structure with small tail at higher binding energy attributed to different types of oxygen containing functional groups in the GO scaffold. The peak-fitting of C 1s spectrum of GO can be deconvoluted into five chemically-shifted components (Figure 3a).

The C 1s peak centered at 284.9 eV as a major component of GO scaffold is attributed to C-C/C-H carbons.

The peak at relatively lower binding energy (284.0 eV) is due to C=C carbon and suggests the presence of unoxidized domains of sp 2 carbon. Three additional overlapped components at 286.8, 288.6 and 289.5 eV are assigned to the hydroxyl/epoxide (C-OH/C-O-C), carbonyl (C=O) and carboxyl (COOH) functional groups, respectively [44][45]. Figure 3 indicates that all characteristic peaks of C 1s spectra remained intact, when Mo 6 clusters were immobilized on the GO nanosheets and no significant changes were observed. Furthermore, the loading of Mo 6 clusters onto GO nanosheets was deduced by appearance of additional peaks at 724.5, 229.3 and 70.2 eV attributed to Cs 3d, Mo 3d and Br 3d, respectively (Figure 4). The Mo 3d peaks at 229.3 and 232.4

eV are due to the split components of 3d 5/2 and 3d 3/2 , respectively and revealed the Mo(II) oxidation state [24,46]. The positions of Mo 3d 4a). The Cs 3d 3/2 peak couldn't appear explicitly, as it was diminished by nearby strong peak of O (KVV) as shown in Figure S1 (Electronic supporting information).

The internal core of the octahedral Mo 6 clusters consists of eight Mo-Br i bonds and the remaining six Br -are placed on apical position through Mo-Br a bonds having strong ionic character. These apical Br -functionalities are prone to interact with other chemical functionalities and facilitate their immobilization onto solid supports.

The XPS spectra of Mo The thermal behaviour of GO, Mo 6 clusters and their composites were examined by thermogravimetric analysis (TGA). GO showed two characteristic major weight losses in the range of 50-110 °C and 170-240 °C, owing to evaporation of trapped water molecules and thermal decomposition of labile oxygen functionalities in GO scaffold, respectively (Figure 5). The Cs 2 Mo 6 Br i 8 Br a 6 cluster is found to be thermally stable upto 500 °C, watt white cold LED flood light. During the photoreduction, 0.5 mL reaction mixture was withdrawn gradually. Then 1 µL withdrawn reaction mixture was injected in a GC-FID system for quantitative estimation.

The performance of the catalysts was evaluated in term of methanol (MeOH) yield as it was the only liquid reduction product. The MeOH yield (µmol.g -1 cat.) with respect to time using Cs cluster units presented in 0.1g amount of catalyst were found to be 0.015 and 0.011, respectively. The TON of methanol per mole of graphene oxide was calculated by considering the equivalent formula of GO as C 6 (O), which was found to be 0.04. The formation of methanol was further confirmed by HPLC analysis as depicted in Figures S5 andS6. For the quantitative estimation of methanol, a calibration curve was plotted by injecting an exact amount of methanol (1 µL) in GC-FID (Figure S8). Then the concentration of methanol produced by photo-induced CO 2 reduction was calculated by injecting 1.0 µL of the final reaction solution in GC-FID and comparing the peak area with calibration curve for methanol standard. Further HPLC analysis was performed to confirm the formation of methanol as the only liquid product obtained from the photoreduction of CO 2 .

Gaseous reaction products were analyzed by GC-FID and GC-TCD and confirmed that there was no gaseous product such as H 2 , CO, CH 4 formed during the photoreduction of CO 2 . As methanol was obtained as the only product from CO 2 photoreduction under the described experimental conditions, so the catalytic selectivity (CS) (conduction band) and sp 3 carbon domains (valence band) on GO sheet. Due to presence of several conductive sp 2 carbon domains and non conductive sp 3 carbon containing domains on GO sheet, a band gap is created and it works like a semiconductor [51][52][53]. The band gap for GO was found to be 2.9 eV in good agreement with the literature values of 2.9-3.7 eV as reported by Hsu et al. [50]. On the other hand, Mo 6 clusters due to their lower band gap values have good visible light absorbance; hence to improve the performance of GO in visible light, Mo 6 clusters were immobilized on the surface of GO via replacement of labile bromine atoms. After absorption of visible light, Mo 6 cluster was transformed to excited state Mo 6 * by HOMO to LUMO transition and the excited Mo 6 * transfers electron to the conduction band of GO (Fig. 9). These electrons in the conduction band of GO were used for the reduction of CO 2 to give methanol [54][55][56]. The positively charged Mo 6 + cluster gets electron from valence band of GO and transforms back into its original (Mo 6 ) state. Water splitting at the valance band of GO possibly provides necessary electrons and protons for the reduction of CO 2 to methanol. Furthermore, due to the better charge mobility and higher surface area, GO improves the catalytic performance synergistically [57][58][59].

6Mo

Fig. 1 :

 1 Fig. 1: A schematic illustration of (a) GO nanosheet decorated with various oxygen functionalities, (b) immobilization of

whereas (TBA) 2

 2 Fig (TB flow

  Fig equ Cs 2

  Fig

Fig

  Fig

  inner and apical ligands and HRTEM image of GO loaded with ample Cs 2 Mo 6 Br i

			8 Br a	6 /
	(TBA) 2 Mo 6 Br i	8 Br a	6 clusters.

Cs 2 Mo 6 Br i 8 Br a 6 / (TBA) 2 Mo 6 Br i 8 Br a 6 clusters on the GO nanosheets, and (c) molecular structure of Mo 6 cluster 5 representing position of

Cs 3d Intensity, cps Binding Energy, eV GO-Cs 2 Mo 6 Br i 8 Br a x GO-(TBA) 2 Mo 6 Br i 8 Br a x 240 236 232 228 224 800 1000 1200 1400 (b) GO-Cs 2 Mo 6 Br i 8 Br a x GO-(TBA) 2 Mo 6 Br i 8 Br a x Mo 3d 3d 3/2 3d 5/2

  e. GO-Cs 2 Mo 6 Br i In order to probe the chemical states and plausible interaction between the Mo 6 clusters and GO, all samples were thoroughly characterized by XPS.

	GO	Fit components Sum Raw spectra		(Cs) 2 Mo 6 Br i 8 Br a 6 GO GO-	(a)	i Br	Sum Raw spectra GO-Cs 2 Mo 6 Br i 8 Br a x (c ii )	(a)	(b)
	at 1466 and 1383 cm -1 associated with C-H bending modes revealed the presence of butyl chains in the (TBA) 2 Mo 6 Br i 8 Br a 6 clusters . The appearance of C-H vibrations in the GO-(TBA) 2 Mo 6 Br i 8 Br x a composite (inset of Figure 2b) suggested the loading of (TBA) 2 Mo 6 Br i 8 Br a 6 on the GO nanosheets. Both GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x composites exhibited the characteristic vibrations of oxygen functionalities of GO, indicating that most of oxygen functionalities were remained intact and presented in both GO-Cs 2 Mo 6 Br i 8 Br a x H 2 O Molecules 294 292 290 288 286 284 282 280 Binding Energy, eV 294 292 290 288 286 284 282 280 80 78 76 74 72 70 68 66 64 62 60 Binding Energy, eV Br 3d % Transmittance, a.u. GO-(Cs) 2 Mo 6 Br i a Fit components 8 Br x O=C-OH C=O C-OH & C-O-C C-C C=C Intensity, a.u. O=C-OH C=O C-OH & C-O-C C-C C=C Br 3d Intensity, a.u. (c i ) Br a 80 78 76 74 72 70 68 66 64 62 60 Intensity, a.u. Br i
	and GO-(TBA) 2 Mo 6 Br i	8 Br a x composites. 3500 3000 GO-(TBA) 2 Mo 6 Br 2500 i 8 Br a x Binding Energy, eV 2000	1500	(c)	1000
									Wavenumbers, cm Raw spectra -1 Sum
									Fit components
									Intensity, a.u.	GO (TBA) 2 Mo 6 Br i GO-(TBA) 2 Mo 6 Br 8 Br a 6 i 8 Br a C=O O=C-OH C-OH & x C-O-C	C-C / C-H	C=C	(b)
	3000 Fig 3: High resolution C1s XPS spectra of (a) GO, (b) GO-Cs 2 Mo 6 Br i 2900 2800 292 290 288 286 284 Binding Energy, eV % Transmittance, a.u. 8 Br a	282 x and (c) GO-(TBA) 2 Mo 6 Br i 280	8 Br a	6 samples. The
	C-H Bending modes deconvoluted components illustrate the various types of oxygen functionalities in the materials. C-H Stretches
									3500	3000	2500	2000	1500	1000
									Wavenumbers, cm -1 (a)
	6500
	8 Br a x and GO-x composite, and (b) GO, 6 composite. The presence of C-H vibrational modes in GO-x clusters and GO-Cs 2 Mo 6 Br i 8 Br a 8 Br a 8 Br a 6 clusters and GO-(TBA) 2 Mo 6 Br i 8 Br a Fig. 2: FTIR spectra of (a) GO, Cs 2 Mo 6 Br i (TBA) 2 Mo 6 Br i (TBA) 2 Mo 6 Br i 8 Br a 6000 Intensity, cps
	(TBA) 2 Mo 6 Br i	8 Br a	6 composite reveals the loading of (TBA) 2 Mo 6 Br i	8 Br a	6 clusters on the GO nanosheets.
				730	728	726	724	722	720	718
									Binding Energy, eV
	the presence of ample hydroxyl, epoxy, carboxylic, carbonyl, phenolic, etc. functional groups in the GO
	scaffold. The Cs 2 Mo 6 Br i	8 Br a x has no vibrational signature, except a strong band at 1620 cm -1 , attributed to the
	bending mode of adsorbed water molecules in the cluster compound (Figure 2a).
									10

x were thoroughly characterized by FTIR, XPS, TG-DTA, HRTEM, EDX and UV-Vis analyses. The FTIR spectrum of GO (Figure 2) exhibited strong characteristics vibrational bands at 3419, 1732, 1622, 1373, 1263 and 1064 cm -1 , which are attributed to O-H stretch, C=O stretch, H 2 O bending modes and C=C stretch, O-H bending and C-O stretch, respectively. The strong intensities of these vibrations revealed

Cs 2 Mo 6 Br i 8 Br a x Cs 2 Mo 6 Br

			i 8 Br a 6
										(d ii )
										(d i )
	80	78	76	74	72	70	68	66	64	62	60

(TBA)

2 Mo 6 Br i 8 Br a 6 GO-(TBA) 2 Mo 6 Br i 8 Br a x Br 3d Br 3d Br

								Br i
			80	78	76	74	72	70	68	66	64	62	60
								Br i
							Binding Energy, eV
	Fig 4: XPS spectra of (a) Cs 3d; (b) Mo 3d; (c) Br 3d regions for Cs 2 Mo 6 Br i	8 Br a x and GO-Cs 2 Mo 6 Br i	8 Br a	x ; (d) Br 3d
	regions for (TBA) 2 Mo 6 Br i	8 Br a	x and GO-(TBA) 2 Mo 6 Br i	8 Br a	x composites.
	XPS analyses of GO, GO-Cs 2 Mo 6 Br i	8 Br a x and GO-(TBA) 2 Mo 6 Br i	8 Br a x samples were carried out to elucidate

a Intensity, a.u.

  5/2 and 3d 3/2 components in the GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x are comparable to those obtained for their pure Mo 6 clusters. These results suggest that the cluster core Mo 6 Br i 8 remained intact during cluster immobilization on the GO nanosheets. The Cs 3d peaks at 724.5 and 739 eV are attributed to the split components of 3d 5/2 and 3d 3/2 , respectively and indicate the presence of Cs(I) in the GO-Cs 2 Mo 6 Br i 8 Br a x composite (Figure

  6 cluster compounds and their immobilized composites GO-Cs 2 Mo 6 Br i 8 Br a x and GO-(TBA) 2 Mo 6 Br i 8 Br a x are shown in Figure 4c-d. The Br 3d peak in Mo 6 clusters is composed of two components at 70.2 and 68.6 eV, which are attributed to covalently bonded inner Br atoms (Br i ) and apical bromide ions (Br a ) of Mo 6 clusters, respectively [47-48]. The Br 3d spectra of GO-Cs 2 Mo 6 Br i Figure 4c-d) showed that the peak component at 68.6 eV, attributed to Br a ligands, is almost diminished, demonstrating the elimination of Br a ligands. This suggests the reaction of Br a of Mo 6 clusters with hydroxyl / oxygen carrying functional groups of GO nanosheets. As a result, Mo 6 clusters immobilized on the GO nanosheets by covalent linkage and Br a are eliminated during immobilization process.

	cluster and GO nanosheets. The experimental ratio of Br/Mo = 2.6 in (TBA) 2 Mo 6 Br i	8 Br a 6 cluster agrees with
	the theoretical value (2.33). After immobilization of (TBA) 2 Mo 6 Br i	8 Br a 6 cluster, the ratio of Br/Mo decreased
	to 1.6, which is closer to the expected value (1.33) for a Br i /Mo cluster core. Recently, Sutherland et al. have
	further demonstrated that labile terminal nitro groups at apical position in Mo 6 clusters were eliminated during
	their chemical attachment to solid supports [49].
	UV/Vis absorption spectra of Cs 2 Mo 6 Br i	8 Br a 6 and (TBA) 2 Mo 6 Br i	8 Br a 6 in acetonitrile exhibited sharp
	peaks at 242 and 243 nm due to the inter ligand transition and small humps at 386 nm and 380 nm due to the
	metal-to-ligand charge-transfer (MLCT) transitions, respectively (Figure S2a,b). For GO, a characteristic
	absorption band at 227 nm due to π→π* transition of discrete units of sp 2 carbon conjugated network and a
	small hump at 300 nm due to n→π* of carbonyl groups were observed. Immobilization of Mo 6 clusters to GO
	enhanced the overall absorption pattern as shown in Figure S2.
	In order to confirm the visible light absorption, we calculated the band gap of GO, Cs 2 Mo 6 Br i	8 Br a 6 cluster, GO-
	Cs 2 Mo 6 Br i	8 Br a x composite, (TBA) 2 Mo 6 Br i	8 Br a 6 cluster and GO-(TBA) 2 Mo 6 Br i	8 Br a x composite with the help
	of Tauc plot as shown in Figure S3. The band gap value for GO was found to be 2.90 eV was in well
	agreement with the reported literature [50]. For Cs 2 Mo 6 Br i	8 Br a 6 and (TBA) 2 Mo 6 Br i	8 Br a 6 clusters compounds,
	two band gap values were obtained; the 4.07 eV and 4.16 eV were due to the interligand transitions while the
	2.40 eV and 2.47 eV were ascribed to the MLCT transitions, respectively. After immobilization of Mo 6 clusters
	on GO, the band gap values significantly decreased to 0.9 and 1.5 eV for GO-Cs 2 Mo 6 Br i	8 Br a x and GO-
	(TBA) 2 Mo 6 Br i	8 Br a x composite, respectively.
					8 Br a x and GO-
	(TBA) 2 Mo 6 Br i x composites (The atomic % of Mo and Br components in both (TBA) 2 Mo 6 Br i 8 Br a	8 Br a 6 cluster and GO-(TBA) 2 Mo 6 Br i	8 Br a X
	composites were calculated based on XPS study, to determine the interaction between the (TBA) 2 Mo 6 Br i	8 Br a
					13

  6 cluster + hν → 6Mo 6 * cluster (excited state)

	6Mo 6 * cluster (excited state) + GO → 6e -(GO conduction band) + 6Mo 6	+ cluster
	6e -(GO conduction band) + CO 2 + 6H +	→ CH 3 OH + H 2 O
	3H 2 O → 6H + + 6e -+ 3/2 O 2
	6Mo 6	+ cluster + 6e -	→ 6Mo 6 cluster
	Overall reaction CO 2 + 6H + +6e -	→ CH 3 OH + H 2 O

showed strong characteristic vibrational modes in the range of 3000-2800 cm - 1 , attributed to the C-H stretches of methylene and methyl groups (Figure2b). Furthermore, strong vibrations


Maria Muguerza

Cyril Briand

Nicolas Jozefowiez

Sandra U Ngueveu

Energy-efficient supplying strategies in assembly lines

Keywords: Supplying strategy, assembly lines, energy consumption, mixed integer linear programming

This paper focusses on the problem of supplying the workstations of assembly lines with components during the production process. The components are assumed to be packed in pallets (or boxes). The pallets are shipped from a co-called supermarket and delivered to the workstations using tow trains with wagons. For each time period, depending on the production plan, a workstation has its specific component demand expressed as a fraction of pallets. The problem is to define the sequence of workstations that have to be visited by the tow train at each time period, as well as the number of pallets delivered at each stop. For that specific problem, this paper presents a Mixed Integer Linear Program that aims at minimizing the energy consumption of the supplying strategy. We particularly discuss which are the most significant energy parameters for that problem and how they can be explicitly modelled into the MILP. We also provide computational experiments and discuss the limitation of the approach.

Introduction

Nowadays, in most factories row materials supplying is governed by the Just In Time (JIT) paradigm. It concentrates on improving the efficiency, service and quality of the production processes without any regard for their energy performances. There is a need of reviewing JIT under an energy-awareness perspective, aiming at balancing economic and environmental performances. Focusing on the internal logistic of factories, the Eco-Innovera project [START_REF]EASY project[END_REF] "Energy-Aware feeding Systems" (EASY) aims at favoring energy-aware practices, taking explicitly energy consumptions factors into account at the very heart of the optimization procedures. Motivated by the challenge associated with the preceding problem, the project target is the design of new optimization and simulation models for the feeding systems of assembly processes, intending to determine which the major energy-consumption factors are and how they can be optimized to reduce energy consumption, still keeping a satisfactory economical efficiency.

In general, feeding systems of assembly lines (see Figure ) are composed by a central warehouse (or a supermarket), several workstations organized in sequence in the production line and a fleet of vehicles (tow trains) that are in charged of delivering the component to the workstations. The components are packed in pallet or boxes. The supermarket is a decentralized area of material supplies, which is located next to the assembly line, where the component boxes are stored.

The time is discretized in a set of tours or delivering periods. For each period, a workstation has a component consumption (possibly periodic) expressed i, terms of boxes. At each tour, the tow trains load the boxes which have to be shipped to the assembly line, follow a route, and stop at the appropriate workstations for delivering its boxes. The routes are usually fixed and start and finish at the supermarket. The number of boxes that a tow train can transport in the same tour is limited. The number of boxes available at each workstation should never exceed the storage capacity of the workstation (which is usually low). A supplying strategy defines the number of boxes that the vehicle has to deliver to each workstation at each period.

supermarket, which makes it easier and more practical by avoiding any interruption since all the requests and deliveries are managed at a single place near the assembly line. In addition, the assembly lines are two-sided where both sides of the line (left and right) are used in parallel. A schematic view of the assembly line with a single supermarket is represented in Figure 1. In a world where natural resources are limited, issues related to energy are becoming more and more important. Vehicles in factories travel a significant quantity of kilometers for supplying the workstations, causing effects in economic and energetic expenses. Whether they use electric energy or fossil fuel, their energetic consumption is not negligible and more and more attention has to be paid for exhibiting energyefficient supplying strategies. Several factors that are inherent to the problem have impact on the energy consumption. We take interest in determining the most significant ones.

Many researches on energy consumption in feeding systems' optimization focussed on eco-driving solutions [START_REF] Saboohi | Model for developing an eco-driving strategy of a passenger vehicle based on the least fuel consumption[END_REF] or considered the inventory routing problem [START_REF] Archetti | Formulations for an inventory routing problem[END_REF]. The eco-driving model identifies optimal driving strategies for the purpose of reducing the consumed energy in different traffic congestion, based on vehicle speed and gear ratio decisions. This kind of model, suitable for urban transportation, seems to be inappropriate for our problem where the driving profile cannot be controlled in this way. Our problem turns out to be closer to an inventory routing problem but without the routing part, the routes being fixed. Based on this observation, Fathi et al. [START_REF] Fathi | A Multiobjective Optimization Algorithm to Solve the Part Feeding Problem in Mixed-Model Assembly Lines[END_REF] proposes a mathematical model in order to minimize the energy consumption in an assembly line taking into account the distance travelled by vehicles.

In this paper, we intend to show that minimizing the travelled distance does not necessarily implies the minimization of the energy. We prove that other parameters can significantly influence the energy spending and should be explicitly taken into account inside the optimization procedures. These parameters are taken from the energy mathematical relationships, which are established in the article written by Garcia et al. [START_REF] Garcia | A power demand simulator of electric transportation systems for distribution utilities[END_REF]. The remainder of this paper is organized as follows. An energy consumption analysis is proposed in Section 2.1 and the most influential energy consumption factors are described in Section 2.2. A mixed linear program for energetic optimization is described in Section 2.

Energy modeling

Many factors influence the energy consumption of vehicles, the most significant being the traction force and the rolling resistance. The effect of stopping a vehicle in a workstation, then restarting it, has also an important effect in terms of energy expenses [START_REF] Deleplanque | The minimization of the number of stops problem[END_REF].

Energy consumption analysis

The mathematical relationships between the physical quantities of energy (E), work (W) and power (P) are :

E = P * dt (1)
P = F * v (2)
E = W = F * dist (3)
E refers to the energy consumed expressed in Joule, P the power consumed in Watt, F makes reference to the force that causes a power consumption and dt, v and dist represent the time diferencial, the speed or velocity and the distance, respectively. The forces that have more influence on the power consumed by the vehicle are: the traction force (F t ) and the rolling resistance (F r ), in Newtons (N). The traction force is used to generate motion between an object and a tangential surface, and it depends on the mass (m T ) and the acceleration of the vehicle (a(t)). The rolling resistance is the force resisting the motion when a body rolls on a surface and varies in function of the load (m T ), the rolling coefficient (C r ) and the gravity (g). The parameter (m T ) represents the mass of the vehicle (determined by the kind of vehicle) in addition to the transported load, which varies according to the demand. Therefore, the expresion that expresses the energy is:

F t = m T a(t) (4)
F r = m T gC r (5
E = m T (a + gC r )v(t)dt (6)
Looking at the energy formula ( 6), all the parameters, except the gravity, depend on the vehicle characteristics. For sake of simplicity, the acceleration (constant) and the maximum speed are assumed known. The acceleration and deceleration are considered equal and have opposite directions. Thanks to the literature, the rolling coefficient between a wheel and the floor is also known. Eventually, the only parameters that vary are the mass transported and the working time of the vehicle.

Energy is the integral of the consumed power over the time, thus it can be easily calculated. Figure 2.1 represents the consumed energy between two workstations. Three phases can be differentiated according to the vehicle state. The first phase corresponds to the acceleration phase where a peak of energy is produced, due to the acceleration. The second phase begins when the speed of the vehicle is constant and there is no acceleration anymore. Finally, in the deceleration phase, the energy consumption is null. The transported load, the travelled distance (length of a tour and number of tours) and the number of stops ( [START_REF] Deleplanque | The minimization of the number of stops problem[END_REF]) change the energy consumption. We need to analyze the influence of each parameter so as to find optimal strategies.

Significant energy-spending parameters

The transported load is directly proportional to the consumed energy [START_REF] Deleplanque | The minimization of the number of stops problem[END_REF], so if the load increases, the consumed energy too. Depending on the way of delivering the load, the energy consumption is different. The load can be supplied in a constant, decreasing or increasing way, whether the box weights are identical, the heaviest boxes are delivered first or, at the reverse, last . The most suitable strategy for the purpose of energy savings is the decreasing one, due to the fact that most part of the distance is travelled with less load. The effect of stopping at a workstation and restarting cause an important boost in the energy demand. During the acceleration phase, a "peak" in the power needed is produced, so a rise in the energy consumption. Decreasing the number of vehicle stops in every workstation reduces the number of acceleration phases, hence the peaks of energy.

The travelled distance is also directly linked to the energy. Nevertheless, doing some experimentations, we exhibit some situations where it is better to travel over a long distance (make more tours) and transport less load, rather than reduce the travelled distance and carrying a high load.

Energy-aware mathematical modelling

In this section, a mixed integer linear programming (MILP) model is presented. This model includes the previous influential factors. The delivered load and the travelled distance are taken into account in the energy formula. Depending on the vehicle stop policy, four scenarios can be distinguished for computing the energy spent between two contiguous workstations. The energy consumption depends whether the vehicle is stopped in the workstation i or not and whether it is going to stop in workstation i+1 or not. Figure 6:

(M i+1 + M v )(E b (d v + 2d a ))
Where M i+1 is the load that is delivered to workstation i+1 and M v , the weight of the vehicle. E b and E a are the energy consumed when the tow train has a constant speed (E b ) and when it is in the acceleration phase (E a ), respectively. The distance travelled during the starting phase denoted d a and the one of the constant speed phase is referred as d v .

This four possibilities can be expressed as a single expression using decision variables. To indicate that the vehicle stops or not in a workstation, a binary variable X i is introduced. It equals 1 only if the vehicle stops in workstation i. The energy consumption between two workstations can be expressed by equation [START_REF] Vehicle Used | [END_REF].

(M i+1 + M v )(E b d v + X i E a d a +(1 -X i )E b d a + (1 -X i+1 )E b d a ) (7)
Extending the energy consumption between two workstations to the whole route and considering all the time periods, the total energy consumed can be computed. It corresponds to the objective function of the following MILP. Each tour has a maximum of n workstations and the total number of possible tours is denoted by N T . Y t indicates whether the tour t is done or not, P it , is an integer variable that represents the number of delivered components in workstation i during tour t and Z ijt is a linearization variable which is enforced, thanks to linear constraints, to be equal to the product X it X jt .

The mixed integer linear programming model is:

Min Z = N T t=1 n i=1 [Y t M v E b (d v + 2d a ) + M v d a (X it (E a -E b ) -E b X i+1t ) n j=i+1 M j (P jt E b (d v + 2d a ) + Z ijt d a (E a -E b ) -Z i+1jt E b d a ) (8)
St:

P it + IL it-1 -d it = IL it ∀ (i, t) (9)
n i=1 P it ≤ A ∀ (i, t) (10)
P it + IL it-1 ≤ C i ∀ (i, t) (11)
X it ≤ Y t ∀ (i, t) (12)
IL it ≥ 0 ∀ (i, t) (13)
P it ≥ 0 ∀ (i, t) (14)
X it ∈ [0, 1] ∀ (i, t) (15)
Y t ∈ [0, 1] ∀ (i, t) (16)
The objective function (8) aims at minimizing the total energy consumption, which is proportional to the number of delivered components P it , the fact of stopping in workstation i during tour t (X it ) and the realization of a tour Y t . The first set of constraints (9) ensures the demand satisfaction during each period and for each workstation. The second set of equations (10) ensures the capacity limit of the vehicle and the third set of constraints (11) the capacity of the workstations. Next set of equations ( 12) ensures the production of a tour, thus if the vehicle stops in a workstation during the tour, Y t is equal to 1. Equations ( 13) -( 16) defines the domain of each variable.

Experimentation

We use the problem instances proposed by Fathi et al. [START_REF] Fathi | A Multiobjective Optimization Algorithm to Solve the Part Feeding Problem in Mixed-Model Assembly Lines[END_REF] where only a single vehicle is considered. For filling our energy model with the vehicle characteristics, we use the parameters from a real tow train [START_REF] Vehicle Used | [END_REF]. The instances vary according to the total number of workstations and tours considered, as well as the component demand profiles and the storage and transport capacities. Since Fathi et al. only pays attention to the travelled distance minimization, the weight of each box and the distance between workstations were not specified. We propose to analyze the impact of those parameters with respect to the energy consumption.

The MILP model has been solved using the mathematical solver Gurobi, which is freely available for academics.The solution gives an energy-optimal supplying plan that specifies the number of components delivered to each workstation at each tour.

We highlight that our model allow to minimize either the total travelled distance or the energy consumption, by simply replacing one objective function by the other. That is the reason why we decided to minimize first the travelled distance then, providing to Gurobi this first initial feasible solution, to optimize the energy consumption. Our purpose is to analyse whether it is worth, with respect to energy spending, to take all the factors into account (delivered load, number of stops and travelled distance) or whether , on the contrary, minimizing the travelled distance (i.e., the number of tours) is enough.

Concerning the computational effort needed for solving the MILP, it has to noticed that, while minimizing the total number of tours is rather easy, tackling the energy objective function is much more time expensive. Gurobi was not indeed able to solve to optimality most of the problem instances in a reasonable amount of time. So, we decided to limit the optimization effort to one hour, then to compare the best solution obtained after one hour with the initial one, aiming at measuring the energy savings. Consequently, the energy gaps that we observe are lower bounds of the best possible savings.

The energy consumption in function of the number of workstations involved in a tour is represented in Figures 789101112. in these Figures, the consumed energy for every different situation is compared.

We consider various scenarios of weights for the component boxes. In the scenario mA the weight of each component box is the same, in scenarios mB and mC the weights were established randomly, allowing higher variations in scenario mC than in mB. Even if box weights are random, we impose that the total demanded weight for each tour is identical in each scenario, which allow to compare the scenarios in energy terms. Concerning the distance between the workstation, we also consider three scenarios dist1, dist2 and dist3. In scenario dist1 the inter-distance. It varies randomly in the other scenarios, allowing higher variation in dist3. In any scenario, the length of a tour is enforced to be constant.

All the experiments show that the consumed energy is always lower when all the parameters are taken into account. Increasing the number of workstations make even the energy savings greater, especially when the distances between workstation vary a lot. The mass of the boxes seem also to have an impact, though for the experiments made so far, the effect of mass variations is less obvious.

Conclusion

We can conclude that taking the transported load, the number of stops and the total travelled distance simultaneously into account is worthy. We propose a MILP formulation that takes all theses parameters into account inside the optimization procedure. Nevertheless, the computational effort required for solving efficiently the corresponding model becomes really huge. Additional researches are needed in order to boost the optimization procedure using either more compact MILP formulation or more advanced optimization mechanisms such as valid inequalities generation, variable fixing techniques, or decomposition approaches (column generation).

Figure 1 .

 1 Figure 1. A straight assembly line view with the supermarket-concept

Figure 1 :

 1 Figure 1: General structure of an assembly line.

  3. In Section 3, experiments are provided and a comparison of our model with the distance minimization model of Fathi et al. is made. Conclusions are finally given in Section 4.

) 14 e

 14 Colloque National AIP-Priméca 2/6 La plagne (73) -31 Mars au 02 Avril 2015

Figure 2 :

 2 Figure 2: Energy profile between 2 workstations.

Figures 3 ,

 3 [START_REF] Fathi | A Multiobjective Optimization Algorithm to Solve the Part Feeding Problem in Mixed-Model Assembly Lines[END_REF][START_REF] Garcia | A power demand simulator of electric transportation systems for distribution utilities[END_REF][START_REF] Deleplanque | The minimization of the number of stops problem[END_REF] shows the four possible scenarios, with their corresponding energy formula.

Figure 3 :Figure 4 :

 34 Figure 3: (M i+1 + M v )(E a d a + E b d v )

Figure 5 :

 5 Figure 5: (M i+1 + M v )(E b (d v + d a ))

Figure 7 :Figure 8 :

 78 Figure 7: Consumed energy by the two models considering a mA load.

Figure 9 :

 9 Figure 9: Consumed energy by the two models considering a mC load.

Figure 10 :

 10 Figure 10: Consumed energy by the two models considering a dist1 distance.

Figure 11 :

 11 Figure 11: Consumed energy by the two models considering a dist2 distance.

Figure 12 :

 12 Figure 12: Consumed energy by the two models considering a dist3 distance.

Table 1 .

 1 Notation

					Route
					Workpiece	Rack	Supermarket
					Rack
					Tow train
					for the MILP model
	I		:		set of references
	i	:		index of references ( i I 	)
	T		:		set of tours
	t	:		index of tours (t T)
	i IL	0 :	inventory level of reference i at the beginning of the first tour
	A	:		capacity of tow train (volume)
	C		:	capacity of correspondence station to reference i (volume)
	i
	d		:	demand of reference i on tour t
	it
	n		:		number of references
	NT	:	number of tours
			:		importance coefficient for the first objective
			:		importance coefficient for the second objective
	p		:	amount of reference i transported on tour t
	it

e Colloque National AIP-Priméca 1/6 La plagne (73) -31 Mars au 02 Avril 2015

e Colloque National AIP-Priméca

La plagne (73) -31 Mars au 02 Avril 2015

e Colloque National AIP-Priméca 6/6 La plagne (73) -31 Mars au 02 Avril 2015


Maria Muguerza

Cyril Briand
email: briand@laas.fr

Nicolas Jozefowiez

Sandra Ulrich Ngueveu
email: ngueveu@laas.fr

Victoria Rodríguez
email: vrodriguez@unav.esmatias.urenda.moris@his.se

Matias Urenda Moris

A mass-flow MILP formulation for energy-efficient supplying in assembly lines *

Keywords:

This paper focuses on the problem of supplying the workstations of assembly lines with components during the production process. For that specific problem, this paper presents a Mixed Integer Linear Program (MILP) that aims at minimizing the energy consumption of the supplying strategy. More specifically, in contrast of the usual formulations that only consider component flows, this MILP handles the mass flow that are routed from one workstation to the other.

Introduction

In general, feeding systems of assembly lines are composed by a central warehouse, several workstations organized in sequence and a fleet of vehicles (tow trains) in charge of delivering the components to the workstations. The components are packed in pallet or boxes. The supermarket is a decentralized area of material supplies, located next to the assembly line. For building up a supplying strategy, time is discretized in a set of delivering periods. For each period, a workstation has a component consumption (possibly periodic) expressed in terms of boxes. At each tour, the tow trains load the boxes which have to be shipped to the assembly line, follow a supplying route, and stop at the appropriate workstations for delivering its boxes. The supplying routes are usually fixed and start and finish at the supermarket. The number of boxes that a tow train can transport in the same tour is limited. The number of boxes available at each workstation should never exceed the storage capacity of the workstation (which is usually low). A supplying strategy defines whether a vehicle has to stop at each workstation at each time period and the number of boxes that should be delivered.

In a world where natural resources are limited, issues related to energy efficiency are becoming more and more important. Vehicles in factories travel a significant quantity of kilometers for supplying the workstations, causing effects in economic and energetic expenses. Whether they use electric energy or fossil fuel, their energetic consumption is not negligible and more and more attention has to be paid for exhibiting energy-efficient supplying strategies. Several factors that are inherent to the problem have impact on the energy consumption. We are interested in determining the most significant ones.

In the literature related to the Vehicle Routing Problem (VRP), some researchers take interest in minimizing carbon dioxide emissions. One of the major contribution is due to Bektas and Laporte [START_REF] Bektas | The Pollution-Routing Problem[END_REF] who present the Pollution-Routing Problem (PRP) as an extension of the VRP with Time Windows. The PRP consists of routing a number of vehicles to serve a set of customers within preset time windows, and determining their speed on each route segment, so as to minimize a function comprising emissions and driver costs. The author propose a MILP formulation that allows to optimize both load and speed of the vehicles.The idea of controlling the vehicle velocity on each route segment is fruitful for improving the energy efficiency in the context of long distance transportation problem. However, in a very local transportation context, as the distances travelled during the acceleration phase becomes non-negligible with respect to the one covered at the maximum speed, other parameters can impact the energy consumption.

The problem considered in this paper can be viewed as a particular Inventory Routing Problem (IRP) . We intend to show that minimizing the travelled distance does not necessarily implies the minimization of the energy. We prove that other parameters can significantly influence the energy spending. The remainder of this paper is organized as follows. An energy consumption analysis is proposed in Section 2. A MILP for energy optimization is described in Section 3.

Energy modeling

The forces that have more influence on the power consumed by the vehicle are: the traction force (F t = m T a(t)) and the rolling resistance (F r m T = gC r ), in Newtons (N). The traction force is used to generate motion between an object and a tangential surface, and it depends on the mass (m T ) and the acceleration of the vehicle (a(t)). The rolling resistance is the force resisting the motion when a body rolls on a surface and varies in function of the load (m T ), the rolling coefficient (C r ) and the gravity (g). The parameter (m T ) represents the mass of the vehicle plus the transported load, which varies along the tour. The expresion of the energy consumption is E = m T (a(t) + gC r )v(t)dt.

For sake of simplicity, the acceleration, the deceleration and the maximum speed are assumed known and constant.Thanks to the literature, the rolling coefficient is also known. Regarding the energy consumed between two workstations. Three phases are distinguished according to the vehicle state. The first phase corresponds to the acceleration phase where a peak of energy is produced, due to the acceleration. The second phase begins when the speed of the vehicle reach its maximum value. Finally, in the deceleration phase, the energy consumption is null.

The travelled distance is directly linked to the energy although it is not the only significant parameter. Indeed, the energy consumption is different depending on the way of delivering the load. The stops at the workstations also have effects on the energy demand. Decreasing the number of vehicle stops in every workstation can reduces the number of acceleration phases, hence the energy.

Energy-aware mathematical modeling

In this section, a mixed integer linear programming (MILP) model is presented. This model integrates the previous influential factors, and similarly to the formulation poposed in [START_REF] Archetti | Formulations for an inventory routing problem[END_REF], takes advantage from a basic flow formulation. Nonetheless, instead of taking the flow in terms of number of components into account, the mass of the shipped components is considered. We assume that only one kind of pallet can be delivered to a given workstation, each having a well-known mass. Therefore, once the mass of delivered components known, the number of components can be easily deduced. Reasoning in terms of masses is interesting since the energy spent for bringing a pallet to one location i to another location j is proportional to its mass. Therefore, one can considered directly inside the MILP formulation the energy cost (C ij ), which represents the energy consumption for shipping one mass unit directly from i to j with j > i.

The component mass brought from i to j during period t is noted M t ij . Decision variables Z t i represent the number of components left at workstation i during period t. They can easily be deduced from the values of the M t ij variables. Eventually, inventory flow variables IL t i , deduced from the Z t i values, are also modelled.

Using the above decision variables, the energy minimization MILP can be formulated as follows. The objective function (1) aims at minimizing the energy consumption, which is proportional to the mass M t ij traversing each arc (i,j ) during period t. Constraints (2) model flow conservation together with demand satisfaction. Constraints (3) ensure that the vehicle capacity A is never exceeded and enforce variables Y t to be set to one when a tour is carried out in period t. The set of equations (4) ensures that the inventory level at workstation i never exceeds the workstation storage capacity c i . Constraints (5) enforce the difference 1 mi ( j<i M t j,i -j>i M t i,j ) to be integral. The constraints (6) impose that the mass brought back to the depot equals the vehicle mass (0 and n + 1 being two virtual nodes associated with the depot). Constraints (7) ensure that, whether some components are delivered in workstation i during period t, the vehicle has to stop in this station at that tour. Set of constraints (8) ensure that, whether the vehicle stops in workstation i at time t, there exist an incoming and an outcoming arc selected at workstation i during period t. Constraints (9) ensure that whether there exists a mass flow between two workstations, an arc between these stations has to be selected too. Equations ( 10)-( 12) define the domain of each variable (m max is the maximum load that the vehicle can transport).

Min z = n i,j N T t C ij M t ij (1)
st:

Z t i + IL t-1 i -d t i = IL t i ∀ (i, t) (2)
n i=1 Z t i ≤ AY t ∀ (t) (3)
Z t i + IL t-1 i ≤ c i ∀ (i, t) (4)
Z t i - 1 m i ( j<i M t ji - j>i M t ij ) = 0 ∀ (i, t) (5)
n i=1 M t in+1 = m v Y t ∀ (t) (6)
Z t i ≤ X t i c i ∀ (i, t) (7) j>i φ t ij = j<i φ t ji = X t i ∀ (i, t) (8)
M t ij ≤ m max φ t ij ∀ (i, j, t) (9)
IL t i , M t ij ≥ 0 ∀ (i, j, t) (10)
Z t i ∈ N ∀ (i, t) (11)
φ t ij , X t i , Y t ∈ {0, 1} ∀ (i, j, t) (12)

Conclusion

We can conclude that taking the transported load, the number of stops and the total travelled distance simultaneously into account is worthy. We propose a MILP formulation that integrates these parameters all together inside the optimization procedure. Nevertheless, the first experiments show that the computational effort required for solving efficiently the model is high. Additional researches are needed in order to boost the optimization procedure using either more compact MILP formulations or more advanced optimization mechanisms such as valid inequalities generation, variable fixing techniques, or decomposition approaches.

* This work was supported by the ECO-INNOVERA-1rst call EASY (ANR-12-INOV-0002).