Céline Platteau

Tatiana Cucu

Bruno De Meulenaer
email: bruno.demeulenaer@ugent.be

Bart Devreese

Marc De

Marc De Loose

Isabel Taverniers

Effect of protein glycation in the presence or absence of wheat proteins on detection of soybean proteins by commercial ELISA

Keywords: Soybean allergy, Maillard reaction, ELISA, protein bound carbonyls, protein 14 aggregation, reactive lysine

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

Introduction

Food allergy is an abnormal immunological response due to sensitization to a food 18 component. With its increasing incidence over the last years, it represents an important health problem, especially in industrialized countries, where it has been estimated to affect about 2% 20 of the adult population and up to 8% of children [START_REF] Monaci | Milk allergens, their characteristics and their detection in food: A review[END_REF]). Together with hen's egg, cow's milk, fish, crustacea, peanut, tree nuts and wheat, soybeans belong to the 'big 22 eight', being responsible for 90% of all allergic reactions (Report of the FAO 1995). It also belongs to the list of 14 allergenic products which have to be declared on the label of a food 24 product when present as an ingredient according to Directive 2007/68/EU (European Parliament andCouncil Directive 2007/68/EC 2007). 26 Soybeans (Glycine max) contain about 35% of proteins and at least 21 proteins have been shown to be IgE-reactive. Soybean allergy can occur as a result of primary sensitization, 28 but can also be related to birch pollen allergy through structural similarities between the birch pollen allergen Bet v 1 and the corresponding protein Gly m 4 in soybeans (secondary allergy) 30 [START_REF] Ballmer-Weber | Soy allergy in perspective[END_REF][START_REF] Mittag | Soybean allergy in patients allergic to birch pollen: 436 Clinical investigation and molecular characterization of allergens[END_REF]). The symptoms associated with soybean allergy are mostly mild, such as the Oral Allergy Syndrome (OAS), however, severe reactions 32 have also been described [START_REF] Foucard | A study on severe food reactions in Sweden -is soy protein an underestimated cause of food anaphylaxis?[END_REF][START_REF] Kleine-Tebbe | Severe oral allergy syndrome and anaphylactic reactions caused by a Bet v 1-430 related PR-10 protein in soybean, SAM22[END_REF][START_REF] Mittag | Soybean allergy in patients allergic to birch pollen: 436 Clinical investigation and molecular characterization of allergens[END_REF]). The incidence of soybean allergy has increased during the last years probably because 34 more products are being produced with soybean [START_REF] Wilson | Allergenic proteins in soybean: Processing and reduction of P34 allergenicity[END_REF]. Soybean is used as an ingredient in a broad range of food products, both for its texturizing and nutritional benefits. 36

Due to its diverse functional properties, such as water and fat absorption, emulsification, foaming and gelation, soybean proteins are often included during food processing steps. 38 However, extensive use of soybean proteins in foods might pose a serious threat to the allergic consumers due to unintentional contamination which might take place during food 40 processing. This can affect the safety of food products and lead to the presence of hidden allergens, therefore making avoidance diets increasingly challenging. 42

Various approaches have been applied to reduce the allergenicity of soybean proteins.

Enzymatic hydrolysis through fermentation is often applied to produce hypo-allergenic 44 formulas, but is strongly dependent on the extent of the hydrolysis degree [START_REF] Frias | Immunoreactivity and amino acid content of fermented soybean products[END_REF][START_REF] Tsuji | Fate of major 458 soybean allergen Gly m BD 30K in rice-, barley-and soybean-koji miso (fermented soybean paste) during fermentation[END_REF][START_REF] Yamanishi | Reduction of soybean allergenicity by fermentation with Bacillus nato[END_REF]. Reduction of soybean allergenicity has also been 46

achieved by allergen conjugation with polysaccharides [START_REF] Arita | Effect of chemical and genetic 388 attachment of polysaccharides to proteins on the production of IgG and IgE[END_REF][START_REF] Babiker | Effect of polysaccharide conjugation or transglutaminase treatment on the 392 allergenicity and functional properties of soy protein[END_REF].

Similarly, the impact of protein modifications occurring as a result of the multitude of 48 chemical reactions occurring during food processing has been studied as well [START_REF] Burks | Identification of Soy Protein Allergens in Patients with Atopic-Dermatitis and 404 Positive Soy Challenges -Determination of Change in Allergenicity After Heating Or Enzyme Digestion[END_REF][START_REF] Codina | Neoallergens in heated soybean hull[END_REF][START_REF] Muller | 446 Allergens in raw and roasted hazelnuts (Corylus avellana) and their cross-reactivity to pollen[END_REF]. In contrast to these studies focusing on the 50 impact of soy protein modification on the antigenicity, only limited studies are available with respect to the impact of these changes. In addition, the impact of these changes on the 52 detectability using immunochemically methods, which are typically used by the industry and enforcement agencies [START_REF] Bando | Quantitative analysis of Gly m Bd 28K in soybean products by a sandwich 398 enzyme-linked immunosorbent assay[END_REF][START_REF] Tsuji | 462 Measurement of Gly-M-Bd-30K, A Major Soybean Allergen, in Soybean Products by A Sandwich Enzyme-Linked-Immunosorbent-Assay[END_REF] has rarely been studied in a 54 systematic manner for soybean protein.. The above-mentioned methods are based on antibody -antigen interactions Therefore, immunoassays are subjected to complications due to the 56 matrix effects, which may result in erroneous results. Moreover, the interaction of food allergens with other food components, which are expected to occur during thermal processing, 58 may lead to changes in the allergen structure and might hinder their detection by the immunological methods. Therefore, it is important to gain knowledge on the influence of food 60 processing on the tertiary structure of food allergens and the changes occurring at the amino acid level as a result of chemical reactions taking place during processing. Moreover, 62 information on the further impact of these changes on the allergenic potential and the detectability of the allergens is invaluable. Van [START_REF] De Lagemaat | In vitro 410 glycation and antigenicity of soy proteins[END_REF]. In this study the impact of glycation in the 66 presence or absence of other proteins on the detectability of soybean proteins by commercial kits was investigated. A similar study on hazelnut proteins already showed that the detection 68 was drastically affected as a result of protein glycation (manuscript submitted).

Materials and methods

Reagents and standards 72

Soluble wheat proteins "Meripro 711" were kindly provided by Syral Belgium N.V.

(Aalst, Belgium). Chemicals and standards of analytical grade were obtained from Sigma-74 Aldrich (Bornem, Belgium), VWR (Leuven, Belgium) and Acros Organics (Geel, Belgium).

Solvents were of analytical grade from Chem-Lab (Zedelgem, Belgium), except for hexane 76 which was of technical grade. Gel and standards for SDS-PAGE were from Bio-Rad (Nazareth, Belgium) while the Krypton protein stain was provided by Thermo Scientific 78 (Aalst, Belgium). The commercial ELISA kits used were Veratox® Soy Allergen (Neogen, Lansing, USA), BioKits Soya Allergen Assay (Tepnel, Flintshire, UK) and Soy Residue 80 (ELISA Systems, Queensland, Australia).

Soybean protein extraction 82

Soybeans were provided by Alpro (Wevelgem, Belgium) and Cargill (Mechelen, Belgium). A mixture was made by taking equal amounts of each kind. Soybeans were frozen 84 with liquid nitrogen and ground in a two-step process, first with a blender (Moulinex, France) and secondly with an Ultra Turrax t25 (IKA, Wilmington, US). 86 Soybean protein extracts were obtained according to the protocol described by Yeung et al. with slight modifications [START_REF] Yeung | Enzyme immunoassay for determination of peanut 472 proteins in food products[END_REF]. Briefly, 100 g of ground beans were 88 defatted with four times 1 L technical hexane. After an overnight evaporation of the solvent at ambient temperature, 10 g defatted powder was mixed with 100 ml of extraction buffer and 90 mixed at 4°C over night. Protein extraction buffer consisted of a 1 M NaCl, 1.47 mM KH 2 PO 4 , 2.7 mM KCl and 10 mM Na 2 HPO 4 . Sodium azide 0.2 g L -1 was added to the 92 extraction buffer to avoid microbial growth. Further a mixture of protease inhibitors (2 mM ethylenediaminetetraacetic acid (EDTA), 5 mM sodium diethyldithiocarbamate (DIECA), 0.5 94 mM benzamidine hydrochloride (BAHC) and 0.2 mM phenylmethylsulfonyl fluoride (PMSF)) was added to the extraction buffer in order to increase the stability of the protein 96 isolates during extraction and storage [START_REF] Bjorksten | Extraction and properties of apple allergens[END_REF]. The protein mixtures were further centrifuged for 1 h at 4°C and 20000 g to separate the lees. Subsequently, the clear 98 supernatants were dialyzed (Spectra/Por 3 membrane tubing, MWCO 3.5 kDa) using the same extraction buffer containing only 0.137 M NaCl, for 48 h at 4°C and then filtered through a 100 0.45 µm membrane filter (Puradisc FP 30,Novolab,Belgium). Finally, the protein solutions were lyophilized and stored at -20°C until further use. This extract will further on in this 102 paper be considered as the reference extract.

Maillard reaction in the absence of soluble wheat proteins 104

Before treatment, soybean protein extracts were dialyzed against the 0.1M phosphate buffer pH 7.4 to remove the nitrogen containing additives (Spectra/Por 3 membrane tubing, 106 MWCO 3.5 kDa ). The reaction systems were prepared in the same phosphate buffer pH 7.4 by making a solution containing 1% (w/v) soybean proteins and 6% (w/v) glucose. The 108 reaction systems in closed falcon tubes were incubated in a water bath at 70°C for up to 48 h.

A reaction system without glucose was used as control. At different time periods, samples 110 were retrieved from the water bath and placed in an ice bath to stop the reaction. Further, the modified protein solutions were immediately used for the determination of reactive lysine 112 residues, protein bound carbonyls and free amino groups. Aliquots for SDS-PAGE were frozen at -20°C until further use. 114

Maillard reaction in the presence of soluble wheat proteins

Two types of model systems were considered. One model system contained 0.25% 116 (w/v) of soybean proteins with 0.75% (w/v) of wheat proteins and 6% (w/v) of glucose in the proteins with 0.5% (w/v) wheat proteins and 6% (w/v) glucose. Respective model systems without glucose were prepared as control. Further these systems were treated similarly as the 120 model systems containing no soluble wheat proteins.

Assessment of the Maillard reaction 122

Determination of protein bound carbonyls

Protein bound carbonyls were measured after converting them to the corresponding 2, 124 4-dinitrophenylhydrazones. Briefly, 300 µl protein solution was incubated with 400 µl 2.4dinitrophenylhydrazine (DNPH) 10 mM in 2 M HCl. Further, the proteins were precipitated 126 with 10% trichloroacetic acid (final concentration) and thoroughly washed with 1 ml ethanol/ethyl acetate (1:1) after derivatization; hence any unreacted DNPH or DNPH bound to 128 any other components than protein was removed. Finally the pellet was redissolved in 8 M urea and absorbance was measured at 370 nm using a Bio-Rad Benchmark Plus microplate 130 spectrophotometer (Nazareth, Belgium). The carbonyl content was calculated using a molar absorption coefficient of 22000 M -1 cm -1 [START_REF] Levine | Carbonyl Assays for Determination of Oxidatively Modified Proteins[END_REF]. 132

Loss of reactive lysine residues

The loss of reactive lysine residues was monitored using derivatization with ortho-134 phthaldialdehyde (OPA) in the presence of 2-mercaptoethanol which yields a fluorescent product with a maximum excitation wavelength of 340 nm and emission of 450 nm using a 136 Spectramax Gemini XPS fluorimeter from Molecular Devices (Brussels, Belgium) [START_REF] Ferrer | Evolution 418 of available lysine and furosine contents in milk-based infant formulas throughout the shelf-life storage period[END_REF]. A standard of β-casein was used to prepare a calibration curve and samples were 138 pre-treated and analyzed using a fluorimeter as previously described [START_REF] Morales | Evaluation of heat-induced changes in 442 Spanish commercial milk: Hydroxymethylfurfural and available lysine content[END_REF].

SDS-polyacrylamide gel electrophoresis (SDS-PAGE) 140

Protein cross-linking or fragmentation was analyzed by SDS-PAGE under reducing conditions. A 4% acrylamide stacking gel and a 10% acrylamide resolving gel were used for 142 protein separation. An aliquot of 10 µL sample containing 6 µg of protein was added to 10 µl Laemmli sample buffer containing β-mercaptoethanol and heated at 95°C for 5 min. The 144 samples were further loaded on the gel along with the molecular weight marker. Precision plus protein unstained standard was used as a reference to determine the molecular weight of 146 proteins of interest. Krypton fluorescence staining was used for bands visualization. The gels were further read using the Molecular Imager PharosFX Plus from Bio Rad (Nazareth,148 Belgium) where the molecular weights of the proteins were recorded.

Protein determination 150

Protein determination was performed by determining the nitrogen content according to the Kjeldahl procedure (AOAC Official Method 981.10 1981). A multiplication factor of 5.71 152 was used to convert nitrogen content to soybean protein content, and a factor of 6.25 to convert to soybean protein content for the mixtures of soybean and wheat protein. The non 154 protein nitrogen (NPN) fraction was determined in the supernatant after a previous protein precipitation with 15% TCA (final concentration). 156

ELISA assays

Samples of the above-described model systems were analyzed using different 158 commercial ELISA kits: Veratox® Soy Allergen, BioKits Soya Allergen and Soy Residue.

These samples include the glucose-free model system containing only soybean proteins 160 (system 1), the model system containing soybean proteins with glucose (system 2), the two model systems containing different ratio's of soybean proteins and soluble wheat proteins 162 (0,5% soybean proteins + 0,5% soluble wheat proteins (system 3) and 0,25% soybean proteins containing different ratio's of soybean proteins and soluble wheat proteins, in the addition of 6% glucose ((system 5) resp. (system 6)). Of the different model systems, samples retrieved 166 after 0, 24 and 48 h incubation at 70°C were analyzed. The untreated sample (= 0 hours of thermal treatment) of the glucose-free model system (1) containing only soybean proteins is 168 referred to as the reference sample. For each sample of the different model systems a dilution series was made (0 -50 -100 -500 -1500 -1500 ng/ml soybean protein) in the extraction 170 buffer/dilution solution of the respective kit. For the BioKits Soya Allergen Assay additional dilutions had to be analyzed. The actual concentrations are the concentrations of soybean 172 proteins based on the absolute protein determination using the Kjeldahl procedure and do not include the wheat proteins that were added in some of the model systems. The dilution series 174 were analyzed in duplo in a similar manner as sample extracts according to the respective manufacturers' manuals. Absorbance measurements were made in a microtiter plate reader 176 (Multiskan MCC/340, Titertek, Huntsville, Alabama). Each kit was calibrated using the respective standard protein solutions included in each kit. The absorbance values were fitted 178 to the protein concentration of the standard solutions, according to a four-parameter-logistic dose-response curve corresponding to following equation (2.1). This equation presents a 180 standard way for the expression of immunoassay data, where a is the maximum dose (protein concentration), b is the minimum concentration, c is the E 50 -value (concentration at half-182 maximal saturation), d is the slope of the curve and conc is the concentration of the standard or sample. 184

( ) a d c conc a b Absorbance +                   + - = 1 (2.

Progress of the Maillard reaction

An increase in protein bound carbonyl amounts as result of the incubation in the 194 presence of glucose has been recorded. A 5 fold increase in protein bound carbonyls was recorded after 24 h incubation of 1% soybean with glucose while after 48 h incubation a 10 196 fold increase was observed (Figure 1). For the model system containing no glucose a 2 time increase was observed. Similarly, an increase in protein bound carbonyls was observed when 198 soluble wheat proteins were added to the model systems. After 48 h incubation with glucose a 7 fold increase in protein bound carbonyls was observed in the case of the model system 200 containing 0.5% soybean proteins. For the model system containing 0.25% soybean proteins a 4 fold increase was observed. 202 Further, changes in the reactive lysine residues have been recorded. About 80% decrease of reactive lysine residues has been recorded upon incubation of 1% soybean 204 proteins with glucose (Figure 2). Interestingly, a 23% decrease in reactive lysine residues after 6 h incubation in the absence of glucose was also recorded. A further decrease in 206 reactive lysine residues was however not recorded, while in the presence of glucose, a steady decrease was obtained throughout the whole incubation experiment. Further, about 90% of the 208 reactive lysine residues were modified in the model system containing 0.5% soybean proteins after 48 h of incubation with glucose while for the model system containing 0.25% soybean 210 proteins almost all the reactive lysine residues were modified. were still slightly visible. Aggregates with masses of over 300 kDa were also formed which hindered their migration into the resolving gel. 230

Detection with commercial ELISA kits

The soybean protein concentration of the standards provided by each individual kit is 232 expressed in parts per million (ppm) which equals ng analyte mg -1 of food product. When testing the commercial kits, the objective was to compare the actual soybean protein 234 concentration of the analyzed samples with the measured concentration, derived from the calibration curve obtained with the kit's standards. To this end, the concentration of the kits 236 standards was converted to the same dimension of that from the model system samples. This For the evaluation of the kits performance, dilution series of the soy protein solutions obtained from the six different model systems were applied in the appropriate concentration 244 range of each kit. The obtained absorbance values were converted to a measured soy protein concentration, which was compared with the actual protein concentration presents in each 246 diluted solution (Figure 4).

Biokits Soya Allergen 248

With the standard solutions included in the BioKits Soya Allergen assay, a calibration curve was constructed according to the four-parameter logistic model, ranging from zero to 250 200 ng mL -1 soybean protein. The soybean protein concentration in the reference samples was overestimated on average by 290% in a concentration dependent manner (Figure 4A). 252 Therefore, at the highest concentrations absorbance values exceeding those obtained in the calibration curve were observed. Upon heating at 70°C without glucose for 24 h the signal 254 decreased with 14% (not shown) and an additional 8% when the treatment was extended up to 48 h. When soluble wheat proteins were added to the model system the average 256 overestimation of the concentration of the samples was diminished to 250%, irrespective the soy/wheat protein ratio. However, at the lowest soybean/wheat protein ratio, a sharper decline 258 in the measured soybean protein concentration was observed upon heating; after 48 h the signal dropped on average 38% (0.75% soluble wheat proteins) compared to an average of 260 27% (0.5% soluble wheat proteins, not shown). Addition of glucose resulted in a further decrease of the measured concentration compared to the samples treated without glucose; the 262 signal was declined with 60% after 48 h heat-treatment. 16 decided to dilute the available standards provided with the kit in order to obtain a more suitable calibration curve. The protein concentration of the native reference sample was 290 underestimated for the samples containing 50 and 100 ng mL -1 soybean protein in a concentration depend manner with on average 12% (Figure 4C). Exceeding the calibration 292 range, the reference samples containing 500, 1000 and 1500 ng mL -1 soybean protein could not be considered and the latter two were not included in Figure 4C. Upon heating without 294 glucose for 24 h at 70°C, only about 2% of the proteins could be detected (not shown) with an additional decrease to only 1% when extended to 48 h. When soluble wheat proteins were 296 added to the system, the measured concentration was lower compared to the reference sample (15 -20%); in contrast to the observations described for the Veratox® Soy Allergen, a lower 298 signal was measured when more wheat proteins were added. Also here an almost complete loss of signal was observed after heating in the presence of soluble wheat proteins, with only 300 1% of the soybean proteins being detected after 24 h (not shown) and 0.2% after 48 h.

Similarly, no additional effect of adding glucose during the heat-treatment could be seen. 302

Discussion

In this study, different commercial ELISA kits for the detection of soy proteins and the 304 influence of protein glycation hereon was investigated. Application of the reference soybean protein solution in the various kits resulted in an extreme overestimation of the actual protein 306 concentration in the Veratox Soy Allergen and BioKits Soya assays, while for the Soy Residue assay an underestimation was recorded. In addition, the amount of soybean proteins 308 that was measured was dependent on the concentration in the sample that was applied. Similar observations have been made in a validation study of five commercial ELISA kits for peanut 310 detection [START_REF] Poms | Inter-laboratory validation study of five commercial ELISA test kits for the 450 determination of peanut proteins in biscuits and dark chocolate[END_REF]. A possible explanation for the differences between the actual and the measured soybean protein concentration that were obtained with the different kits 312 could be found in the nature of the standard solutions of the respective kits. In the Soy Residue kit a single protein, namely the Soy Trypsin Inhibitor, is used as the target analyte. 314

Possibly the factor used to convert the soy trypsin inhibitor concentration to total soy protein concentration is not the same as the actual proportion soy trypsin inhibitor present in the 316 soybean protein extract used in this study. Strategically, the choice for a single marker protein for the detection of soybean proteins in processed food products is especially risky in terms of 318 susceptibility to protein modification and its impact on the detection, as was moreover demonstrated in this study. The high overestimation of the actual soy protein with the BioKits 320 Soya Allergen assay demonstrate that the antibodies that are used to develop the assay recognize better the proteins that are present in the soybean protein extract used in this study 322 than the ones present in the standard solutions of the kit. This could be due to the use of different soybean varieties to produce the protein extract. Immunochemical activity of peanut 324 proteins has also been shown to depend on the peanut variety that is used [START_REF] De Meulenaer | 412 Development of an enzyme-linked immunosorbent assay for peanut proteins using chicken immunoglobulins[END_REF]. These observations indicate that the test results heavily depend on the assay that is 326 Upon heating the reference soybean protein solution in absence of glucose, a decreased detection was observed in the BioKits Soya Allergen, however, still overestimating the actual 330 soybean protein concentration. In contrast, both the Veratox Soy Allergen and the Soy Residue assay lacked total robustness. The complete loss of detectability can probably be due 332 to the denaturation of soybean proteins especially because only minor changes in terms of formation of protein carbonyls and losses of lysine residues were observed as compared to the 334 changes observed in the model systems containing glucose (Figure 1&2). However, the observed loss of intensity of the Gly m 5 allergen upon heating without glucose, indicating 336 denaturation, might explain the lowered detection by the BioKits Soya Allergen Assay [START_REF] Kitabatake | Denaturation Temperature of Soy Protein Under Low Moisture Conditions[END_REF]). The Gly m 6 allergen, with a denaturation temperature of about 90°C 338 could then account for the remaining detection. The Soy Trypsin Inhibitor shows a thermal transition at 60-65°C which could explain the decreased detection with the Soy Residue assay 340 after heating at 70°C [START_REF] Roychaudhuri | Stability of the allergenic soybean Kunitz trypsin inhibitor[END_REF]). Although the denaturation of this protein has been shown to be reversible, fast cooling of the model systems, by putting the samples on ice 342 after heat-treatment, could have prevented a proper refolding.

When soybean proteins are present in a food product, protein cross-linking can occur 344 with other (bulk) proteins, which could have an additional impact on the detection. Therefore the impact of heating soybean proteins in the presence of soluble wheat proteins was 346 investigated as well. In a previous study on hazelnut proteins it was seen that soluble wheat proteins can have an additional impact (manuscript submitted). However, no additional 348 impact on the detection of soybean proteins was seen when soluble wheat proteins were added to the model systems. 350 During heating of the reference soybean proteins with glucose either with our without soluble wheat proteins, changes on molecular level took place as evidenced by a number of 352 parameters investigated (Figure 123). As a consequence, these changes on protein level led to further decreases in detectability especially for the BioKits Soya Allergen Assay. 354 Nevertheless, in the same kit, additional protein cross-linking with other soluble wheat proteins via the produced protein bound carbonyls during glycation, did not further affect the 356 detectability of the soybean proteins. As a complete lack of robustness of the Veratox Soy Allergen and the Soy Residue was observed upon heating soybean proteins alone, no 358 additional impact of their heating with glucose either in the presence or absence of soluble wheat proteins could be seen. 360

The results presented show that the detection of soybean proteins by means of commercial ELISA kits is strongly influenced by processing. This implies that certain assays 362 could falsely declare a processed food product as "free of soybean proteins", while soybean proteins might still be present. Consequently, this holds a risk for the allergic consumer. In the 364 same time some other kits, such as BioKits Soya Allergen and Veratox Soy Allergen, show strong overestimation especially for the native soybean proteins which presents risks for the 366 food industry where the product withdrawal from the market is based on the results obtained with similar commercial kits. In conclusion, there is still a lack of sensitive analytical methods 368 for the detection of traces of allergens in processed foods and thus controls cannot be carried out in a reliable way. Moreover, in this study, by inducing Maillard reaction in soybean 370 proteins in a buffered system, abstraction of possible consequences of protein insolubility and eventual extractability was made. This implies that, within a real food product, it might be 372 expected that estimations of the contamination levels would be even more severely affected. It should therefore be emphasized that the tested ELISA assays can be useful to give an 374 indication of the presence of contaminants, rather than an exact quantitative value of the values ±SD of three independent determinations. □ 1% soybean proteins with 6% glucose; ○ 0.5% soybean with 0.5% soluble wheat proteins and 6% glucose; ∆ 0.25% soybean and 492 0.75% soluble wheat proteins and 6% glucose; ◊ 1% soybean proteins Figure 2 Loss of available lysine due to incubation of soybean proteins in the presence or in 494 the absence of soluble wheat proteins with glucose, data points represent mean values ±SD of four independent determinations. □ 1% soybean proteins with 6% glucose; ○ 0.5% soybean 496 with 0.5% soluble wheat proteins and 6% glucose; ∆ 0.25% soybean and 0.75% soluble wheat proteins and 6% glucose; ◊ 1% soybean proteins 498 Figure 3 SDS-PAGE pattern of the soybean proteins incubated in the presence or absence of wheat proteins with glucose Lanes 1-3: 1% soybean proteins non-treated, 24 h and 48h 500 incubated resp.; lanes 4-6: 1 % soybean proteins and 6 % glucose non-treated, 24 h and 48 h incubated resp.; lanes 7-9: 0.5 % soybean proteins with 0.5 % wheat proteins and 6 % glucose 502 non-treated, 24 h and 48 h incubated resp.; lanes 10-12: 0.25 % soybean proteins with 0.75 % wheat proteins and 6 % glucose non-treated, 24 h and 48 h incubated resp.; MWM -504 molecular weight marker

  phosphate buffer pH 7.4. Another model system contained 0.5% (w/v) soybean 118

  75% soluble wheat proteins (system 4), respectively) and the two model systems 3 and 4 164

  obtained within the range of the respective standard curve were 186 considered in the results. In case of the Veratox® Soy Allergen kit, the kit's procedure prescribes to fit the data to log-logit function, which is a linearization of the four-parameter 188 logistic function. For the Soy Residue and the BioKits Soya Allergen Assay no suggestion is made as to which function to use to fit the data. Data above or below the calibration range are 190 indicated by an asterisk (*) or a 'zero' (0), respectively, in the charts.

Figure 3

 3 Figure3shows the SDS-PAGE pattern of the native and glycated soybean proteins. 212

  was done according to formula 3.1, considering the weight of the sample (mg) extracted by a 238 particular volume of extraction buffer (mL)

  the need for reference materials for the development and validation of standardized methods. 328

Figure 1

 1 Figure captions 488

Figure 4 Figure 1 Figure 2 Figure 3 Figure 4 Figure 4 Figure 4

 4123444 Figure 4 Ratio of measured over actual soybean protein concentration in different ELISA kits 506

  de Lagemaat et al. previously reported 64 reductions of the antigenicity of soybean proteins after in vitro glycation using an in-house

	developed ELISA assay
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	5

http://mc.manuscriptcentral.com/tfac Email: fac@tandf.co.ukFood Additives and Contaminants

Acknowledgements

This work was funded by the Belgian Science Policy (BELSPO) (contract n° 380 SD/AF/03A). The authors would like to thank the Neogen Corporation for kindly providing us the Veratox Soy Allergen ELISA kits for this study. 382

Veratox® Soy Allergen 264

A calibration curve was established with the standard solutions provided in the Veratox® Soy Allergen kit. The standard consists of Soy Protein Isolate and covers a range 266 from zero to 4000 ng mL -1 after conversion of the concentration expressed in ng g -1 , according to formula 3.1. As can be observed from Figure 4B, the actual soybean protein 268 concentration present in the reference solution and its dilutions was highly overestimated on average by 300%. Therefore, the obtained absorbance values for the reference samples with 270 1500 and 1000 ng mL -1 soybean protein were beyond the range of the calibration curve and were not included in Figure 4B. Upon the presence of soluble wheat proteins in the model 272 system the overestimation was even increased up to on average 400% (0.5% soybean protein + 0.5% soluble wheat proteins, not shown) or on average 600% (0.25% soybean proteins + 274 0.75% soluble wheat proteins). However, already after heating the reference sample for 24 h in the absence of glucose the signal disappeared completely. In the presence of soluble wheat 276 proteins only 1% (0.5% soluble wheat proteins) or 8% (0.75% soluble wheat proteins) of the soybean proteins was detected after 24 h heating (not shown) and extension to 48 h lead to a 278 complete loss of signal. Heating the proteins in the presence of glucose confirmed the lack of robustness of the kits. 280

Soy Residue

In the Soy Residue assay the soybean trypsin inhibitor is used as the analytical target. 282

The kit includes standard solutions containing 5, 2.5, 1 and 0 ppm soybean protein corresponding to 500, 250, 100 and 0 ng mL -1 soybean protein after conversion which were 284 used to construct a four-parameter logistic calibration curve. The 5 and 2.5 ppm standards returned equal absorbance values, amounting 2.8-3, consequently quantification can only be 286 performed up to a level of 2.5 ppm. In addition, the 1 ppm standard, being the lowest standard containing the analyte, returned a quite high absorbance value of 2.5. Therefore, it was 288


Bernat Gacias

Julien Cegarra

Pierre Lopez

Scheduler-oriented algorithms to improve human-machine cooperation in transportation scheduling support systems

Keywords: Decision support system, vehicle routing problem, work domain analysis, control modes, user-oriented algorithms, model inversion

A decision support system designed to enhance human-machine interaction in transportation scheduling is proposed. We aim to integrate human factors and ergonomics from the beginning of the design phase and to propose a system fitted with enough flexibility to be able to deal with the characteristics of a dynamic context such as transportation scheduling. In this interdisciplinary approach, a link is done between problem solving methods (operations research techniques and data classification algorithms) and human-machine interaction (solving control modes). A set of scheduler-oriented algorithms favouring human-machine cooperation for problem solving is proposed. Some of these algorithms have been efficiently tested on instances of the literature. Finally, an original framework aiming to assist scheduler in constraint relaxation when the problem becomes infeasible is proposed and evaluated.

Introduction

Nowadays, many firms emphasize the need to support routing performance because of an increasingly competitive environment. For this reason, routing problems have been a predominant application area for decision support systems ever since they were first introduced, as evidenced by bibliometric analyses [START_REF] Eom | DSS Applications Development Research: Leading Institutions and Most Frequent Contributors (1971[END_REF][START_REF] Eom | A Survey of Decision Support System Applications (1988-1994)[END_REF][START_REF] Eom | Survey of Decision Support System Applications (1995-2001)[END_REF]. The starting point for these support systems is usually to tackle the hard combinatorial problem arising from the numerous constraints that have to be taken into account, then to allocate any remaining components to the human planner. Sanderson noted back in 1989 that these support tools are not especially dominant in practice, and very little has changed since then. Vehicle routing is often performed by a single experienced individual or, more rarely, by a small group of individuals (see [START_REF] Cegarra | A cognitive typology of scheduling situations: A contribution to laboratory and field studies[END_REF]). In field studies, humans appear to be crucial for taking the large set of constraints into account and adapting to changes in the domain [START_REF] Sanderson | The human planning and scheduling role in advanced manufacturing systems: an emerging human factors domain[END_REF][START_REF] Jackson | A new model of scheduling in manufacturing: Tasks, roles, and monitoring[END_REF]. At the same time, the integration of transportation technologies (GPS, EDI, GIS) is gradually changing the nature of these decision-making processes, in that vehicle routing systems now combine human planners and technologies.

In this paper, we begin by discussing constraint processing and the importance of modelling the work domain in order to assess planners' constraint processing. After that, the work domain analysis of the vehicle routing problem and the decision support system architecture is introduced at the end of Section 2.

Typical support systems focus on supporting the "mechanical" process of generating the schedule. Algorithms therefore function like a black box: the human has to provide the necessary information and has to adapt the outcome if the route contains errors, but during the actual generation process, s/he plays no role whatsoever. In Section 3 we present how to take into account the cooperative aspects between human and computer, which are both participating in the decision-making process.

Support needs above all to consider how human perform tasks and their performance. In the last section a model inversion framework allowing to support the human in relaxing constraints is presented.

Background

There have been a number of attempts to define the vehicle routing task structure on the basis of hierarchical task analysis [START_REF] Rahimi | A hierarchical task model for dispatching in computer-assisted demand-responsive paratransit operation[END_REF] and cognitive task analysis [START_REF] Wong | Analysing ambulance dispatcher decision making: Trialing emergent themes analysis[END_REF]. However, when [START_REF] Cegarra | A comparison of task analysis methods for planning and scheduling[END_REF] compared the amount of information produced by these two methods, as well as by work domain analysis, they found that the latter was far more exhaustive in identifying constraints, not least because it provides a generic view of constraints and does not focus on usual or known tasks.

Work domain analysis (WDA) was developed by Rasmussen and colleagues [START_REF] Rasmussen | Cognitive Systems Engineering[END_REF][START_REF] Vicente | Cognitive Work Analysis: Toward safe, productive, and healthy computer-based work[END_REF]. Instead of focusing directly on the tasks being considered by the decisionmaker, it looks at the constraints imposed on behaviour by the environment. For an interface designer, this sometimes requires a change of point of view, as Vicente (2000) (p.63) pointed out: "A task can be defined as the set of actions that can or should be performed by one or more actors to achieve a particular goal. In contrast a work domain is the system being controlled, independent of any particular worker, automation, event, task, goal, or interface". As previously stressed, for the purposes of assessing and supporting decision-making, completeness in identifying domain constraints is highly desirable in a vehicle routing system.

WDA is usually performed using an abstraction hierarchy which depicts the constraint space. The higher levels of the abstraction hierarchy describe functional information about the domain, whereas the lower levels describe physical information. Moreover, WDA usually looks at five levels of abstraction [START_REF] Naikar | Work domain analysis: Theoretical concepts and methodology[END_REF]: functional purposes, abstract functions, generalized functions, physical functions, and physical forms. In addition to this breakdown into physical and functional aspects, a part-whole distinction can also be made, taking several levels of details into account: system, subsystem and components. However, this part-whole distinction is not always made, as we will see below.

Initially, WDA was applied to "causal" systems guided by physical laws, as in nuclear power stations [START_REF] Itoh | An ecological interface for supervisory control of BWR nuclear power plants[END_REF], conventional power stations [START_REF] Burns | Navigation strategies with ecological displays[END_REF] and cement milling plants [START_REF] Van Paassen | New visualisation techniques for industrial process control[END_REF]. In "causal" systems, the objective reality is imposed on the human operator [START_REF] Vicente | Cognitive Work Analysis: Toward safe, productive, and healthy computer-based work[END_REF], as opposed to "intentional" systems, where the operator is the main agent of the domain and there are fewer references to an external environment. This is the case of routing problems, in which it is difficult to enumerate the domain constraints because they result from conventions, organizational objectives, formal or informal rules and operators' goal. [START_REF] Vicente | Cognitive Work Analysis: Toward safe, productive, and healthy computer-based work[END_REF] stressed that WDA should be performed independently of usual or known tasks, the aim being to provide an exhaustive breakdown that is resistant to changes in the situation. To this end, we sought to enhance the identification of domain constraints by extending the scope of the domain. More specifically, instead of focusing on currently known or usual cases in one particular situation (which inevitably leads to the inclusion of more details about the current situation), we set out to identify constraints from different situations documented in the literature. Each variant of the VRP would allow us to increase the completeness of the constraint space by con-sidering constraints that might potentially help to identify planners' degrees of freedom. Our search for VRPs in scientific databases yielded more than ten thousand articles. However, only a few of them suggested genuinely new variants of the generic VRP (see [START_REF] Toth | The vehicle routing problem[END_REF]). Instead, researchers tended to develop algorithms to solve known variants or design algorithms for multiple variants (e.g., [START_REF] Pisinger | A general heuristic for vehicle routing problems[END_REF]).

In the main, the variants extend the domain by providing constraints related to a temporal perspective, including customers' time windows (VRPTW) and vehicles with limited capacity (CVRP), and more complex variants, such as trucks and trailers (TTRP), pick-up and delivery (VRPPD), or ones related to multiple depots (MDVRP), which are summarized in Table 1. As previously indicated, the WDA abstraction hierarchy organizes the planners' problem space in terms of different concepts that planners can then use for reasoning within a work system. We also consider that the classic five levels of the abstraction hierarchy are needed to describe the constraint space. Figure 1 summarizes the VRP domain according to these different levels. [START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF]. The decomposition remains very generic according to studies and problem descriptions found in the literature. Nevertheless, most of the vehicle routing problems are covered and only some particular variants may required other components not considered on our analysis (see [START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF] for a scope analysis). A support system architecture aiming to cover some of the major issues concerning humanmachine cooperation is proposed. An important clue often pointed by researchers is the advisability of sharing a unique reference system between the human and the machine. To this purpose, the work domain analysis has been considered crucial in order to identify the information needed to solve the problem. In our support system, as shown in Figure 2, the human and the solving mechanism manipulate both the same entities (physical objects and technical constraints) identified in the analysis. The system is composed by a set of human-machine interfaces seeking to efficiently assist the scheduler to solve problem tasks (see Gacias et al. (2010b) and [START_REF] Cegarra | Implications of technological changes in vehicle routing systems for planners' constraint processing[END_REF] for a detailed description of these interfaces). The information displayed has been identified in the work domain analysis. The idea is to connect these interfaces and the humansystem interaction tools to the user-oriented algorithms presented in the next sections in order to favour human-machine cooperation for problem solving.

Variant Name

Description

Objectives, constraints under consideration Real-life example VRP Vehicle routing problem

A set of vehicles is to be routed from a single depot to multiple customers.

The objective is to minimize the number of vehicles and total travel distance. Workload balancing sometimes appears as a secondary objective of the problem.

See more specific cases below.

Customers-related variants DVRP Dynamic VRP

A VRP where a subset of parameters of the problem (new customers, etc.) evolve in real time.

The objective is to minimize deviations from the original schedule (stability criterion) and minimize total travel distance and the number of vehicles.

Travelling repairman, emergency services, taxicab services. (e.g., [START_REF] Rahimi | A hierarchical task model for dispatching in computer-assisted demand-responsive paratransit operation[END_REF]) VRPTW VRP with time windows A time window is associated with each customer (the interval at the depot is named the scheduling horizon).

The objective is to minimize the size of the vehicle fleet and of the total travel distance and waiting time needed to supply all customers at the required times.

Waste collection (e.g., [START_REF] Kim | Waste collection vehicle routing problem with time windows[END_REF]).

VRPPD VRP with pick-up and delivery

The demands of the customers are indifferently pickups or deliveries. Some variants of the problem are:

• Simultaneous PD: pick-up and delivery for each customer must be carried out simultaneously.

• Mixed PD: customers either have a pick-up or a delivery demand.

The objective is to minimize the size of the vehicle fleet and total travel distance, with the restriction that the vehicle must have sufficient capacity for transporting the commodities to be delivered and those picked up from customers and taken back to the depot. One of the most widely used strategies for solving the problem is delivery first and pick-up second.

Beverage industry

(where filled bottles are delivered while the empty ones are collected), on-demand transportation (in reality, customer time windows are also considered).

ARP

Arc routing problem

The customers' demands are on arcs.

The objective is to minimize the number of vehicles and total travel distance.

Trash collection (e.g., Santos et al.

(

)) Vehicle-related variants CVRP Capacitated VRP 2008
The vehicles are identical and based at a single depot, and only the capacity restrictions for the vehicles are imposed.

The objective is to minimize the number of vehicles and total travel distance, with the restriction that vehicle capacity must be respected.

Milk collection from supplier farms (e.g., [START_REF] Basnet | FleetManager: a microcomputer-based decision support system for vehicle routing[END_REF]).

DVRP

Distanceconstrained VRP

The length (or duration or cost or number of customers) of each vehicle's route is bounded by a prescribed amount.

The objective is to minimize the size of the vehicle fleet and total travel distance.

The travelling auditor (e.g., [START_REF] Mendoza | An evolutionary-based decision support system for vehicle routing: The case of a public utility[END_REF])

HVRP

Heterogeneous fleet VRP

There exist a number of heterogeneous vehicle types (capacity, fixed cost, variable travel cost).

The objective is to minimize the total cost of the routes. The best vehicle fleet composition has to be determined.

Feed compound delivery (e.g., [START_REF] Ruiz | A Decision Support System for a Real Vehicle Routing Problem[END_REF]).

TTRP Truck and trailer RP

The fleet is made up of trucks and trailers. Some customers can only be served by a single truck but others can be served either by a single truck or by a truck pulling a trailer.

The objective is to minimize the total distance travelled, or the total cost incurred by the fleet. The truck's uncoupling and re-coupling of its trailer is authorized in some locations.

Depot-related variants MDVRP Multiple depot VRP

If the customers are clustered around depots, then the distribution problem should be modelled as a set of independent VRPs. However, if the customers and the depots are intermingled, then a multi-depot VRP should be solved.

The objective is to minimize the size of the vehicle fleet and total travel time, and the total demand of commodities must be met from several depots.

Gas and oil station delivery, ready-mixed concrete distribution (e.g., [START_REF] Matsatsinis | Towards a decision support system for the ready concrete distribution system: A case of a Greek company[END_REF]).

Decision support system algorithms

The algorithms integrated in the solving mechanism are presented in this section. The section starts with a formal description of the vehicle routing problem and of the problem constraints. The control modes for problem solving are then presented. Finally, we describe the proposed algorithms.

Formal problem statement

The vehicle routing problem has already been formally defined (see for example [START_REF] Toth | The vehicle routing problem[END_REF] for a collection of linear programs proposed in the literature for different variants of the problem). Our problem is the vehicle routing problem taking into account the constraints identified during the work domain analysis (see Section 2).

A set {C i } i=1..nc of customers has to be served by a set {V j } j=1..nv of vehicles. A vehicle starts its route from a depot in {De t } t=1..nde and ends at the same depot. To each vehicle is assigned a driver in {D l } l=1..ndr . Each customer demands an amount d k i for each product in {P k } k=1..np and a time service ts i .

We describe here some characteristics of problem constraints. First, the customer demands can be indifferently deliveries (d k i > 0) or pick-ups (d k i < 0). Second, each vehicle j is capacity constrained: a maximal weight capacity (C w j ), a maximal volume capacity (C v j ) and a maximal authorized length (C l j ). As the same way, each product k is defined by a weight (P w k ), a volume (P v k ) and a length (P l k ). Let us define the variable x j i to indicate if a customer i is served by the vehicle j.

x j i = 1 if customer i is served by vehicle j 0 otherwise
Two different sets of constraints can be then defined, one for the deliveries the other one for the pick-ups (depending on the sign of the demand), as follows:

nc i np k x j i × |d k i | × P w k ≤ C w j ∀ j = 1..nv nc i np k x j i × |d k i | × P v k ≤ C v j ∀ j = 1..nv (1) max i (x j i × P l k ) ≤ C l j ∀ j = 1..nv
For each vehicle/route a limit for travelled distance (V D j ), time (V T j ), and for the number of customers (V C j ) may also be defined. Let us define R j as the route of vehicle j, D(R j ) and T (R j ) are the travelled distance and the duration of R j , respectively. The following new set of constraints is then defined: [r t , d t ]) and driver (T W Dr l = [r l , d l ]) time windows are also considered. The vehicle departures and arrivals have to take place inside depot time windows. Drivers cannot work outside their time windows. Precedence and immediate precedence constraints between customers can be also defined. If customer i precedes customer i (i ≺ i ), customer i has to be served before customer i , and if customer i immediately precedes customer i , customer i has to be served immediately after customer i by the same vehicle. Finally, allocation constraints between physical objects (depot-vehicle, customer vehicle, product-vehicle, etc.) are also considered.

D(R j ) ≤ V D j ∀ j = 1..nv T (R j ) ≤ V T j ∀ j = 1..nv (2) nc i x j i ≤ V C j ∀ j = 1..nv
Two different interfaces are proposed for problem modelling (see Figure 3). These interfaces are designed from the physical levels (purpose-related functions, physical functions, and physical objects) of the abstraction hierarchy.

User-oriented algorithms: control modes

The role of the human in problem solving has traditionally been a conflict issue between operations research and human factors researchers. Indeed, the latest claimed the advisability of give scheduler a more important role on problem solving process. Function allocation to humans and algorithms in planning and scheduling has been recently raised in van Wezel et al. (2010).

The authors propose principles for function allocation based on the analysis of problem subtasks, human capabilities, and the characteristics of the available algorithms to solve the problem.

In that context, a solving mechanism with different control modes to perform the problem tasks is proposed. Problem solving has been divided into three subtasks: vehicle selection, customer allocation, and route selection. This task division is justified in [START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF] after the interview of two different companies transportation planners. In van Wezel et al. (2010), the authors propose five control modes: manual, advisory, interactive, supervisory, and automatic. The control mode specifies the degree of user participation in problem solving process. In the manual control mode, all decisions are made by the human (none algorithmic assistance is provided). In the advisory control mode, human makes the decisions and an algorithm checks the decisions feasibility. In the interacting control mode, the decision-making process is shared between the human and the algorithms. In the supervisory control, the algorithm is first executed, then necessarily informs the user, deciding to accept or reject the algorithm decisions. Finally, in the automatic control mode, the decisions are completely made by the algorithm (the user is then completely out of the decision-making process).

For problem solving, we think that problem tasks require both scheduler and algorithms participation. Indeed, the complex nature of problem constraints, some of them cannot even be considered on problem modelling phase, demands scheduler participation in order to ensure their satisfaction. On the other hand, problem complexity and the huge amount of computing requirements demands the articulation of transportation technologies (GPS, EDI, GIS) with efficient algorithms in order to obtain satisfactory solutions in a reasonable time. In a real-life context, the full-manual control mode is unreasonable in front of too many items to manage. Indeed, it involves exceeding the cognitive capacities of humans, especially for decision-making under stress. Moreover, a dynamic environment implies painful re-computing. On the other hand, the consideration of the automatic control mode falls out of the scope of our paper focusing on user-oriented algorithms. In that context, we propose a three-phase solving mechanism with different control modes (advisory, interactive, and supervisory) facilitating scheduler participation in decision-making process and integrating efficient algorithms for problem solving and constraint satisfaction.

Vehicle selection algorithms

The goal of vehicle selection algorithms is to assist scheduler to select the vehicles serving customers. Let us define y j as the variable to define if a vehicle j is used to serve customers.

y j = 1 if
d k i × P p k ≤ nv j=1 y j × C p j (3) nc i=1 np k=1 d k i × P v k ≤ nv j=1 y j × C v j (4) max i,d k i = 0 (P l k ) ≤ max yj >0 (C l j ) (5)
If capacity and allocation constraints are satisfied and the number of vehicles proposed by the user is greater than LB nv , the solution is then accepted. If it is not the case, a constraint relaxation assistance based on model inversion techniques is proposed to the scheduler (see Section 4).

Supervisory and interactive control mode: vehicle number minimization

Both control modes focus on vehicle number minimization. In that context, several solutions satisfying capacity and allocation constraints and minimizing the number of vehicles are proposed to the scheduler. Besides, the lower bound LB nv has to be also respected by these solutions.

The problem can be defined as a linear program as follows:

min nv j=1 y j (6) subject to nc i=1 np k=1 d k i × P w k ≤ nv j=1 y j × C w j (7) nc i=1 np k=1 d k i × P v k ≤ nv j=1 y j × C v j (8) nv j y j ≥ LB nv (9) y j ∈ {0, 1} (10)
An exact and polynomial algorithm to solve the problem when either one of the capacity constraints (weight or volume) are relaxed is proposed. Algorithm 1 illustrates volume constraint relaxation. First, a solution satisfying the weight constraint and minimizing the number of vehicles is proposed. Then, the satisfaction of the lower bound is checked. The algorithm starts selecting the vehicle with the biggest capacity among the available vehicles except for the last vehicle (when the sum of vehicle capacities is enough to satisfy the constraint) where the vehicle with the smallest capacity satisfying the constraint is selected.

No polynomial algorithm exists when all constraints (weight and volume constraints) have to be considered. Actually, the problem is known to be NP-hard [START_REF] Garey | Computers and intractability: A guide to the theory of NP-Completeness[END_REF]. In that case, the idea is to solve both relaxations separately in order to find a solution satisfying all capacity constraints. Note that when no feasible solution is found after solving both relaxations (on volume and weight constraints) using Algorithm 1, the integer linear problem (6-10) has then to be solved.

Finally, a list of interesting solutions with a minimal number of vehicles is proposed to the scheduler. This set of feasible solutions are found by trying to replace in the solution each Step 2: Select the vehicle V s with a maximal weight capacity if C w Vs < weight then

Step 3: Add vehicle V s to the solution else

Step 4: Select vehicle V s with a minimal weight capacity but greater than weight

Step 5: Add vehicle V s to the solution Step 6: Update the weight, weight ← weight -C w Vs else

Not feasible problem if number of vehicles of the solution is less than BI nv then

Step 7: Select a vehicle V s

Step 8: Add vehicle V s to the solution vehicle of the original solution by a new vehicle with smaller capacity. Of course, the solution after replacing vehicles has to satisfy the capacity constraints of the problem. The list of solutions proposed is not a complete enumeration; still, a list of interesting solutions is provided to the planner.

The interactive control mode follows the same principle. The vehicle number minimization begins instead from a partial solution proposed by the scheduler. The algorithms propose a set of solutions taking into account the decisions already made by the scheduler.

Again, in case of infeasibility because not enough vehicles are available to satisfy the problem constraints, a constraint relaxation assistance is then proposed (see Section 4).

Customer allocation algorithms

The goal here is to determine for each customer of the problem an allocation on a vehicle. The system offers the possibility to manually allocate the customers to vehicles (advisory control mode); in this case the feasibility of each decision is checked by an algorithm. An algorithm for customer allocation is also integrated in solving mechanism. The goal of the algorithm is to propose a complete solution for supervisory control mode and to complete a partial solution proposed by the scheduler for interactive control mode.

Advisory control mode: scheduler decision checking

An algorithm checking allocation feasibility is proposed. The idea is to verify after each allocation decision that it exists a feasible solution. The algorithm is divided in two parts: a constraint satisfaction checking and a feasible solution searching. A feasibility test using energetic reasoning [START_REF] Lopez | Consistency enforcing in scheduling: A general formulation based on energetic reasoning[END_REF] principles is proposed for the first part. In a scheduling context, Gacias et al. (2010a) have proposed an algorithm integrating setup times in the feasibility test. Indeed, travel times between customers in vehicle routing problem can be seen as sequence-dependent setup times. A limited discrepancy search (LDS) algorithm [START_REF] Harvey | Limited discrepancy search[END_REF] is used to find a feasible solution. If no solution is found, the scheduler has to backtrack to her/his previous choices or has to modify her/his last decision. Further details of the algorithm are described in [START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF].

Supervisory and interactive control mode: customer allocation algorithm

In the supervisory control mode, a complete solution is proposed by the system. The scheduler can however modify the decisions of the algorithm. The sweep algorithm [START_REF] Gillet | A heuristic algorithm for the vehicle dispatch problem[END_REF] principle is used for customer allocation. First, customers that have to be allocated to a specific vehicle because of allocation constraints are allocated. Then, each customer (selected from the sweep algorithm order) is allocated on the non-empty vehicle with the smallest mean distance between the customer and the already allocated customers. If there is not a not-empty available vehicle to allocate a customer, then the customer is allocated to a new vehicle.

The interactive control mode uses the same algorithm to complete a solution partially proposed by the scheduler. The scheduler specifies here the allocation for some of the customers and the algorithm ends the customer allocation from the partial solution.

Route selection algorithms

The goal of route selection is to determine for each route a sequence of customers. Besides the control modes proposed here to solve the problem, the system also integrates several algorithms for solution optimization.

Figure 4 shows the proposed user interfaces for route selection.

Advisory control mode: scheduler decision checking

The scheduler can manually propose a sequence for each vehicle through the human-machine interaction tools. Once again, an algorithm tests the feasibility of the scheduler decisions. First, different dominance rules are tested to verify that the proposed sequence is not inconsistent. After that, if the sequence is correct, a tree search algorithm based on limited discrepancy search is launched to find a feasible sequence with all the customers allocated to the vehicle respecting the subsequence proposed by the scheduler. If no solution is found, the scheduler has to backtrack to her/his previous choices or has to modify her/his last decision.

Supervisory and interactive control mode: customer allocation algorithm

The idea here is to propose a feasible complete solution for the problem. To this purpose, different algorithms are proposed: the first algorithm is based on classical customer insertion techniques traditionally used in the literature and the other two algorithms are proposed for solution optimization. The goal is to propose first a feasible solution with the first algorithm (a solution with a structure respecting schedulers' criteria). The scheduler may then use the optimization algorithms to try to improve the quality of the solution. The use of metaheuristics and constraint programming has already been proposed efficiently to solve the vehicle routing problem [START_REF] Pesant | A constraint programming framework for local search methods[END_REF][START_REF] Caseau | A meta-heuristic factory for vehicle routing problems (meta-programming for meta-heuristics)[END_REF][START_REF] Backer | Solving Vehicle Routing Problems using Constraint Programming and Metaheuristics[END_REF]. Once customer allocation is performed, the algorithm uses the principle of savings and regrets algorithms proposed by [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF] and [START_REF] Liu | A route-neighborhood-based metaheuristic for vehicle routing problem with time windows[END_REF], respectively. Each customer is inserted in the position of the sequence that minimizes route length. If a customer cannot be inserted in a sequence, a tree search based on LDS is launched in order to find a feasible sequence. The drawback of the savings algorithm is to find a feasible solution when the problem is over-constrained [START_REF] Kilby | A comparison of traditional and constraint-based heuristic methods on vehicle routing problems with side constraints[END_REF]. It is worth remarking the importance of having an efficient insertion algorithm, thus each time a customer cannot be inserted, a feasible sequence has to be found and the tree search may be time consuming. We then propose to integrate in the algorithm the regrets criterion favouring first the sequence of the most conflicting customers. The proposed algorithm then sequences the customer following the savings criterion except when it exists a customer with only two or less feasible positions into the sequence; in that case the conflicting customer is sequenced first.

Two solutions optimization algorithms are also proposed. The goal of one of them is to independently optimize the routes of the solution. The second algorithm is a local search algorithm for complete solution optimization. The execution time for each algorithm is specified by the scheduler.

The first algorithm is a Climbing Discrepancy Search (local search method using a limited discrepancy search for neighbourhood exploration) [START_REF] Milano | On the relation between complete and incomplete search: an informal discussion[END_REF]. The goal of the algorithm is to optimize the customer sequence for one route. Customer allocation cannot be therefore modified during the optimization process. The main advantage of the approach is that problem solution may be improved by optimizing only some routes and keeping the same structure for the remaining part. Indeed, in the process of solving problem the scheduler may define a solution structure according to her/his criteria; it is then advisable to propose optimization algorithms disturbing as less as possible the solution structure.

The second algorithm for problem solution optimization is described in Algorithm 2. The idea is to integrate the algorithm for route optimization in a classical local search scheme. For each iteration, a customer reallocation based on geographical criteria is first proposed and then the modified routes are optimized using the CDS algorithm. We define now the reallocation operators used in the algorithm. We find some interesting operators in the literature such as the Relocation and the Exchange proposed by [START_REF] Savelsbergh | The vehicle routing problem with time windows: Minimizing route duration[END_REF] or the CROSS-exchange defined by [START_REF] Taillard | A tabu search heuristic for the vehicle routing problem with soft time windows[END_REF]. Our goal is to find as soon as possible the best customer allocations. The idea is to focus on geographical criteria to define our neighbourhoods. The first operator uses the distance between the customers and the gravity centre of the routes to define the new allocations. The distance between a customer and the gravity centre of its route is compared with the distance between the customer and the centroid of the rest of routes. The idea is to determine the customers (C i ) nearer of the centroid of another route (G Rj ) than its own route centroid (G Rq ). If it is the case, the customer C i is then selected to be allocated to route R j . The second operator compares for each customer the distance between the customer and the other customers of the problem with the distance between the customer and its own route centroid. If it exists a customer C i of route R q nearest of the customer C i of route R j than its own route centroid G Rq , the customer C i is then selected to be allocated to route R j . Each new allocation can be defined as {C i , R q , R j }, where C i specifies the customer that will be reallocated from route R q to route R j . We propose then to explore the next neighbourhoods:

• (1) Each new reallocation. The neighbourhood consists in allocating customer C i on the vehicle of the route R j .

• (2) Each allocation exchange between two customers. If it exists an allocation defined as {C i , R q , R j } and another allocation with this form {C i , R j , R q }, then allocations of both customers are modified at the same time. The customer C i is allocated to the vehicle of route R j and the customer C i is allocated to the vehicle of route R q .

It is worth to note our neighbourhoods are sub-neighbourhoods of Relocation and Exchange operators proposed by [START_REF] Savelsbergh | The vehicle routing problem with time windows: Minimizing route duration[END_REF].

Algorithm 2: Solution optimization algorithm

Step 1: An initial solution Sol is proposed by scheduler or solving mechanism algorithms foreach route of Sol do

Step 2: Optimize routes independently using the CDS algorithm (authorised discrepancies are limited) if a better solution BetterSol is found then

Step 3: Update the solution, Sol ← BetterSol while execution time defined by the scheduler is not achieved do

Step 4: Determine all reallocations {C i , R q , R j } defined for each operator foreach new reallocation defined by neighbourhood (1) and(2) do

Step 5: Optimize route R q using the CDS algorithm

Step 6: Optimize route R j using the CDS algorithm if a better solution BetterSol is found then

Step 7: Update solution and go back to Step 4, Sol ← BetterSol

Computational results

The efficiency of the solution optimization algorithm (Algorithm 2) is evaluated here. The algorithm was coded in C++ and run on a 2.8 GHz personal computer with 3.8 Go of RAM under the Linux Debian operating system (Version 5.0.7). The algorithm was tested for the well-known instances of the literature (25, 50, and 100 customers) proposed by [START_REF] Solomon | Vehicle routing and scheduling with time window constraints: Models and algorithms[END_REF] for the VRPTW.

Table 2 shows the results of the comparison between our algorithm (called LS+CDS) and the best-known solutions for distance minimization. The best-known solutions are the optimal solutions for most of the small instances (25 and 50 customers). Computation time is limited to 120 seconds for small instances (25 and 50 customers) and it is increased to 300 seconds for the instances with 100 customers. The first column NV specifies the average number of vehicles, DIST defines the average travelled distance, MeanDev is the mean deviation to the best-known solutions, NBest is the number of instances where our algorithm reaches the bestknown solution and finally TimeBest indicates the average time to find the best solution. These results show that the algorithm is very efficient for small-size instances. Indeed, in most cases good or even optimal solutions are reached very quickly for the instances with 25 and 50 customers. Besides, the mean deviation with the best-known solutions is also acceptable for all instances. It is clear that the performance of the algorithm decreases with instances size. However, the goal of the algorithm is to propose a good performance to optimize a given solution (sometimes defined by the scheduler). Our algorithm is then hardly dependent of the initial solution contrarily to the most efficient metaheuristics of the literature. Anyway, it is always interesting for robustness to keep as much as possible the structure of the initial solution, specially if the scheduler participates in the design of the solution.

Best

Model inversion: A new approach for constraint relaxation

One of the weak aspects of decision support systems is the lack of mechanisms proposed to deal with infeasible problems. Obviously, an infeasible problem is unsolvable, that means it does not exist a solution satisfying all the problem constraints. It is then necessary to relax some constraints in order to find a solution. It is important to stress the advisability of involving the scheduler in constraints relaxation process. Indeed, in planning and scheduling context, the scheduler has a wide knowledge of problem constraints and, more importantly, s/he also has the ability to negotiate constraints relaxation with other problem actors (customers, drivers, etc.).

In the literature few approaches have been proposed to deal with infeasible problems. On the one hand, human-machine interaction experts propose to focus the research on the way the problem constraints are presented in order to facilitate constraint relaxation process to the scheduler [START_REF] Higgins | Interaction in hybrid intelligent scheduling[END_REF](Higgins, , 2001)). On the other hand, [START_REF] Jussien | e-constraints: explanation-based constraint programming[END_REF] proposes, in a fully automated context, algorithmic mechanisms based on the concept of "explanation" to determine the constraints to relax to get a feasible problem.

In the next section, an original framework is proposed aiming to support scheduler by offering a list of interesting possibilities for constraint relaxation making the problem feasible.

Model inversion principles

The study is focused on how to manage constraint relaxation when the problem is not satisfiable. The idea is to propose model inversion mechanisms in order to identify the constraints to relax and how these constraints have to be modified. The model inversion consists in the exchange of roles between the decision variables and the problem fixed parameters. In model inversion, decision variables become parameters restricting decision space and parameters become decision variables that can be used to deduce inferences. The main idea is then to determine how problem parameters have to be relaxed in order to get a feasible problem.

Figure 5 shows the proposed user interface for vehicle selection implementing model inversion mechanisms.

Figure 5: Model inversion user interface for vehicle selection. The scheduler selects the family of parameters authorised to be modified for constraint relaxation when the problem is infeasible (in that case, the weight-capacity constraint is violated).

The first step for model inversion development is the identification of the parameters integrating the constraints. Once the parameters have been identified, it is then necessary to propose a model inversion mechanism adapted to each parameter for each problem constraint. These mechanisms are launched after a constraint is violated (see Figure 5). First, it is necessary to determine the parameters susceptible to be relaxed. Indeed, there always exist parameters having a more important influence on the violated constraint. In a second phase, an algorithm computes how the parameters have to be relaxed in order to satisfy the violated constraint.

We stress that the complicated part of model inversion framework are the design of efficient algorithms to deal with a set of parameters. For instance, it is sometimes difficult to select between a set of customers which one has to be postponed for another day or to decide the customer which its forecasted service time has to be modified in order to get a feasible solution.

In that context, model inversion techniques based on data classification methods in order to determine the most suitable constraints to relax in priority are proposed. However, these generic algorithms need to be adapted to each problem constraint to increase their efficiency.

Data classification based algorithms

In this section, the generic data classification methods with different constraint adapted criteria used to propose to the scheduler a list of options for constraint relaxation are described. Data analysis offers a set of methods designed to structure data information in order to identify connections between individuals. The goal of data classification methods is to propose sets of homogeneous individuals. A similarity measure is often used to group individuals. The goal of the measure is to quantify the similarity between two individuals. The groups of individuals are then proposed trying to maximize the similarity between the individuals.

Two data classification methods to group the customers following a geographical and a temporal criterion are considered. When a constraint is violated, an analysis of the groups of clients in order to determine the most suitable constraints to relax is proposed. However, we remark that these methods can be easily extended to other relevant properties of individuals, in our case the customers, in order to adapt the algorithms to each problem constraint.

Geographical criterion

We propose the k-means algorithm [START_REF] Forgy | Cluster analysis of multivariate data: efficiency versus interpretability of classifications[END_REF][START_REF] Lloyd | Least square quantization in PCM[END_REF] to group customers following a geographical criterion. The k-means algorithm is used to group the nc customers into K clusters (P = P 1 ∪ P 2 ∪ . . . ∪ P K ) as homogeneous as possible in relation to a similarity measure defined for each couple of customers. In the context of a geographical criterion, the customer location is selected as similarity criterion. The properties defining customers are then their location (x i , y i ). The number of clusters is the number of vehicles of the solution. Thus each cluster may be seen as a route.

The similarity measure is then the distance between the customers:

d ij = (x i -x j ) 2 + (y i -y j ) 2 (11)
The k-means algorithm is an iterative algorithm (see Algorithm 3). First, the K means are created. The initial means are uniformly spread along the longest axis of the problem taking the centroid of the customers as the centre point. Then, at each iteration the distance between the customers and the means is computed. Each customer is allocated to the cluster with the nearest mean. Finally, the K means are updated as the centroid of the new clusters. The algorithm stops when no more changes between the clusters of two successive iterations are observed.

Algorithm 3: k-means algorithm for geographical customer classification

Step 1: Select a set of K means repeat

Step 2: Compute distances between each customer and the means (Equation 11)

Step 3: Allocate each customer to the cluster with the nearest mean

Step 4: Update the means of each cluster P k until no changes are observed between the clusters of two successive iterations;

Once the sets of customers are performed, an analysis of the clusters in order to determine the best constraints to relax is performed. This analysis has however to be adapted to the constraint being violated.

Case study: Postpone a customer delivery

We study here the case when a customer delivery has to be postponed because the problem constraints cannot be satisfied. In the geographical criterion context, the customers being part of the same cluster have a high probability to be served by the same vehicle. In that case, a parameter allowing to identify the less homogeneous customers for the clusters is proposed. This parameter (dm i ) is calculated using the mean distance between the customer and the other customers of its cluster and the distance between the customer and the nearest depot in order to penalize the customers not located close of the depots.

dm i = min l∈Depot d il + j∈ P k d ij |P k | (12)
The customer can be then classed according to this parameter. A priority list of customers to postpone based on geographical information can be proposed to the scheduler .

Temporal criterion

A temporal-based customer classification is also proposed. In that case, the clustering approach uses the Dynamic Cluster Algorithm (DCA) introduced by [START_REF] Diday | Une nouvelle méthode en classification automatique et reconnaissance des formes : la méthode des nuées dynamiques[END_REF]. DCA is an extension of k-means algorithm. The similarities here can no longer be shown as an Euclidean distance; an allocation function (a dissimilarity measure) and also a new way to represent the clusters (the centroid does not exist anymore) need to be defined. A dissimilarity measure based on customer time windows [START_REF] Levy | Méthodes par décomposition temporelle et problèmes d'ordonnancement[END_REF] is proposed. The parameter measures the degree of centring between two time windows. Figure 6 displays the main relations between two time windows and the value of the dissimilarity measure (δ) for each configuration. The dissimilarity measure is defined as follows:

1 0 enclosed not-enclosed centered disconnected adjoined i δ(i, j) -min(d j -r i , d i -r j ) d i -r j d i -r j d i -r j d j -r i d j -r i d j -r i i i i i j j j j j
δ(i, j) =              1 - min(dj -ri,di-rj ) max(dj -ri,di-rj ) if min(d i , d j ) ≥ max(r i , r j ), 1 - min(dj -ri,di-rj ) 1 n n i=1 (di-ri)
otherwise.

An allocation and a representation function are also necessary. The core of a cluster P k is the customer C k , which its release date r k is the nearest to the middle of the time interval [min j∈P k r j , max j∈P k r j ]. At each iteration, a customer C i is allocated to the cluster P k with a core C k minimizing δ. The algorithm is very similar to k-means algorithm (Algorithm 3); at each iteration the distances between the customer and the core of the K clusters (δ(C k , C i ), ∀ k = 1..K) are computed for each customer, each customer is then allocated to the cluster with a minimal distance (minimize δ) and finally the core of the new clusters are updated.

To initialize the cores, Algorithm 4 is proposed. The cores are uniformly distributed all along the release times scale. Step 4: Determine the core of the first cluster C 1 , C 1 ← customer i with r i nearest to t 1

Step 5: Initialise cluster counter, k ← 2 for (k ≤ K) do

Step 6: Compute rough time for cluster k, t k ← t k + step

Step 7: Determine the core of cluster k, C k ← customer i with r i nearest to t k and the customer is not already a core of a cluster

Step 8: Increase cluster counter, k ← k + 1

Case study: Postpone a customer delivery The dissimilarity measure groups together customers with time windows interaction. Once the groups of customers are performed, the clusters are analysed in order to identify the most conflicting customers. For each cluster, a parameter called critical index (CI k ) is defined. This parameter is computed as the ratio between the number of customers of the cluster and the number of vehicles. The critical index gives an idea of how critical the cluster is. Indeed, a big value for the CI k means that the planner may have some troubles to serve the customers of the cluster because of customers time windows incompatibilities. The CI k is then used as a parameter to point the customers to postpone in case of infeasibility.

Algorithm 5: Identification of a customer to postpone the delivery for a temporal-based criterion

Step 1: Core initialization for temporal-based classification (Algorithm 4) repeat

Step 2: Compute δ between each customer and the K cores

Step 3: Allocate each customer on a cluster which its core minimizes δ

Step 4: Update cores of each cluster P k until no changes are observed between two successive iterations;

Step 5: Identify clusters of conflicting customers (CI k > limit) for each conflicting cluster do

Step 6: Launch algorithm for geographical classification (Algorithm 3)

Step 7: Save into a list L the customer with the most important dm i Step 8: Propose the customer of the list L with the biggest dm i as the first customer to postpone

We observe that travel time between customers is not considered in temporal-based classification. To cover this lack, both classification criteria (geographical and temporal) are mixed in order to increase efficiency. In that case, the idea is first to target a set of critical customers using the CI k index and then the k-means algorithm is used to select a customer between the critical customers. Algorithm 5 describes the mechanism to propose a priority list of customers to be postponed based on temporal information.

Computational results

The algorithms based on data classification methods proposed to select the best customer to postpone the delivery are evaluated in this section. These algorithms are first compared on a set of small-size instances (9 customers) where a complete enumeration of all feasible solutions can be performed in a reasonable time. The number of feasible solutions is used to evaluate the flexibility provided by the criteria. Indeed, a larger number of feasible solutions widen planners' degrees of freedom facilitating the construction of real-world adapted solutions. The instances of the second set are more realistic instances with 25 customers. Model inversion criteria have been compared considering a geographical optimization function: travel distance minimization.

The small-size instances (9 customers) of the capacitated vehicle routing problem with customer time windows (CVRPTW) are optimally solved. Indeed, each instance is solved nine times suppressing each time one customer from the set of customers. The results obtained taking out of the problem the customer for which the criteria give the priority are compared: the GDC criterion (mean distance between customers of the same cluster after the geographical-based clustering), the GDD criterion (largest distance with the nearest depot), and the TDC criterion (mean distance between customers of the same cluster after the temporal-based clustering).

To generate the instances, customers have to be eliminated of the CVRPTW small-size instances of [START_REF] Solomon | Vehicle routing and scheduling with time window constraints: Models and algorithms[END_REF]. Depending on customer locations, Solomon's instances are classified in three groups: clustered customers (C), random customers (R), and mixed customers (RC).

In Table 3, the first column NbOptDist specifies the number of instances the optimal solution for distance minimization is reached when the problem is solved without the customer proposed by the criterion. The total number of instances is put in brackets. The second column (AvgDev ) specifies the average deviation from the optimal solution for the instances where the optimal solution is reached when the suppressed customer is not the customer selected by the criterion. AvgPos indicates the average position of the solution when the solutions are sorted in a non-decreasing order of the objective function. For example, if the solution reached when the problem is solved without the customer selected by the criterion is the second best solution, then its position is 2. Finally, the last column (NbSol ) represents the number of instances the decision of the criterion is the decision providing a bigger number of feasible solutions.

These results show that GDC is the most efficient criterion for small instances: it reaches the best suitable solution over a half of the instances. The AvgPos is around 3 for the GDC criterion that means that when the criterion decision is not the best suitable decision, the proposition of the criterion keeps being a good proposition. The GDC and GDD criteria outperform the TDC criterion, which is not a surprise because the objective (distance minimization) used to evaluate the algorithms is a geographical-based function such as GDC and GDD, contrarily to TDC which is a temporal-based criterion. The last line shows that all criteria together are really efficient. Indeed, the most suitable choice is selected for 37 over 56 instances by one of the three criteria when distance minimization is considered. It is interesting to note that following the advices of the criteria generally (38 over 56 instances) leads to a problem accepting a greater number of feasible solutions, thus increasing planners' opportunities for their behaviour in order to find better real-world suited solutions. The results are very similar than for the instances with only 9 customers. The GDC criterion proves to be the most efficient criterion, especially for instances with clustered and randomized customers, type C and type R, respectively. Nevertheless, these results show that GDD criterion is the best criterion for the RC instances. It can be due to a possible loss of efficiency of the kmeans algorithm for this type of problems (the clusters of customers after the k-means algorithm do not correspond exactly with the vehicle routes determined by the problem-solving algorithm). However, it is worth remarking that the efficiency of all criteria considered together is really good. Indeed, for most of instances (42 over 56) there is a criterion that proposes a decision reaching one of the five best solutions for distance minimization and the average deviation is less than 5 % meaning that the criteria complement each other.

Conclusions

In this paper, a decision support system for transportation scheduling has been proposed. An interdisciplinary approach has been followed for the system design; a link has been done between problem solving methods (operations research techniques and data classification algorithms) and cognitive engineering techniques (interfaces design, decision sharing between human and machine). The paper focuses on the solving mechanism and scheduler-oriented algorithms integrating the decision support system.

First, a work domain analysis using an abstraction hierarchy for the vehicle routing problem has been described. A new approach is proposed, based on the study of the different variants of the problem detailed in the literature. This approach facilitates to reach completeness in domain constraints identification at the same time that the abstraction hierarchy remains event-and device-independent. The main advantage of the approach is that the problem representation after the analysis does not restrict strategies for scheduler behaviour. The domain constraints revealed by the analysis make up the reference system shared by the human and the machine. In that context, a decision support system architecture has been then defined around the work domain analysis.

The solving mechanism designed to assist the scheduler to perform each problem subtask has been presented. To this purpose, different solving control modes seeking to facilitate scheduler participation in problem solving processes have been proposed. The idea behind control modes, contrarily to a fully automated approach, is that system, through user-interaction tools, gives the schedulers the possibility to built solutions according to their preferences, letting to take into consideration new constraints and/or occasionally violating some of the constraints in order to find better real-world adapted solutions. In that context, a set of useful user-oriented algorithms to assist the scheduler to solve the problem have been proposed and efficiently tested on instances of the literature.

Finally, an original framework aiming at supporting problem constraints relaxation when this becomes infeasible has been proposed. Constraint relaxation support is a topic barely discussed in the literature. In this paper, the principles of a methodology for constraint relaxation seeking to support the scheduler by proposing a list of interesting possibilities for action have been defined. In that context, model inversion mechanisms using data classification algorithms have been proposed and efficiently tested on instances of the literature. The results show the interest of such a Human-Machine cooperation in terms of number of feasible solutions obtained and their quality. Besides, it is worth remarking that the proposed model inversion mechanisms are generic enough to be used on different kinds of problems others than vehicle routing. Indeed, spatial and temporal decomposition methods suit and can be easily adapted for planning and scheduling problems.

Figure 1 :

 1 Figure 1: Work domain analysis of the generic vehicle routing problem Details of the entities corresponding at each abstraction level are presented in[START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF]. The decomposition remains very generic according to studies and problem descriptions found in the literature. Nevertheless, most of the vehicle routing problems are covered and only some particular variants may required other components not considered on our analysis (see Gacias

Figure 2 :

 2 Figure 2: Decision support system architecture

  Then, time windows (T W C i = [r i , d i ]) are considered for each customer. Each customer has to be served inside the interval of its time windows. Besides, depot time windows (T W De t = (a) User interface for objects information (b) User interface for physical objects constraints

Figure 3 :

 3 Figure 3: User interfaces for problem modelling. The scheduler specifies the problem components identified by the work domain analysis: physical objects (vehicles, drivers, goods, depot, and customers), their appearance characteristics (capacity availability and capacity required), their location characteristics (routing network), the temporal constraints (driver and customer time windows and customer service times), and the problem constraints between the objects (allocation and precedence constraints).

Algorithm 1 :

 1 Number of vehicles minimizationStep 1: Compute the weight to transport, weight ← is bigger than weight to transport, nv j C w j ≥ weight, then while (weight > 0) do

Figure 4 :

 4 Figure 4: User interfaces for route selection. One of the interfaces is a geographical representation allowing the scheduler to visualize routes on a map and displaying relevant information about the problem constraints. User interaction tools are also available facilitating the scheduler participation on problem solving.

Figure 6 :

 6 Figure 6: Main configurations between two time windows

Algorithm 4 :

 4 Core initialisation for temporal-based classification Step 1: Compute temporal horizon T H, T H ← max i=1..nc r i -min i=1..nc r i Step 2: Compute step, step ← T H K Step 3: Compute rough time for first core, t 1 ← min i=1..nc r i + step 2

Table 1 :

 1 Variants of the vehicle routing problem

  vehicles required to solve the problem (LB nv ). The lower bound is computed from different problem information. Customer and depot time windows and route constraints such as maximal route time, distance or number of customers are therefore considered to compute LB nv . Further details about lower bound computation are described in[START_REF] Gacias | Une approche interdisciplinaire pour l'ordonnancement des transports[END_REF].

	nc	np
	i=1	k=1

vehicle j is used to serve customers 0 otherwise

Advisory control mode: user solution checking

The advisory control mode consists in a user solution checking. It verifies whether a solution proposed by the user is a feasible solution. To this purpose, global capacity constraints (weight, volume, and length) are first verified (Equations 3, 4, and 5). Then, allocation constraints involving vehicles such as vehicle-customer or vehicle-product are checked in order to ensure the presence of a given class of vehicle necessary to satisfy customer demands. Finally, the number of vehicles proposed by the user to serve customers ( nv 1 y j ) is compared with a lower bound of minimal number of

Table 2 :

 2 Comparison of solution optimization algorithm with the best-known results

			-Known			LS+CDS
	Instances	NV	DIST	NV	DIST	MeanDev NBest TimeBest
	C1-25	3	190.59	3	190.59	0.00 %	9 (9)	0.56
	R1-25	4.92	463.37	5	474.58	3.70 %	4 (12)	0.58
	RC1-25	3.25	350.24	3.25	354.99	9.24 %	7 (8)	0.15
	C1-50	5	361.69	5	361.69	0.00 %	9 (9)	15.38
	R1-50	7.75	766.13	8.5	805.2	5.24 %	0 (12)	15.03
	RC1-50	6.5	730.31	6.87	768.86	7.19 %	3 (8)	2.21
	C1-100	10	826.7	10	873.98	7.36 %	4 (9)	49.56
	R1-100	13.25 1173.61 14.50 1305.32	11.87 %	0 (12)	184.10
	RC1-100 11.12 1341.33 13.37 1410.14	5.39 %	0 (8)	42.67

Table 4

 4

	displays the results for the instances with 25 customers. A complete enumeration

Table 3 :

 3 Results for distance optimization and the number of feasible solutions of solutions in a reasonable time is not possible here, the instances have been therefore solved using Algorithm 2. The last column of Table3calculated from the complete enumeration of feasible solutions has no longer any sense in Table4. For this reason, the column NbSol has been replacing by Nb5Best indicating, in that case, the number of instances the decision proposed by the criterion is among the five best decisions considering distance optimization.

	56 instances
	nc = 25	NbOptDist AvgDev AvgPos Nb5Best
	GDC-C	6 (17)	2.74 %	6.55	9
	GDC-R	4 (23)	4.48 %	6.84	15
	GDC-RC	3 (16)	11.35 %	9.85	6
	TDC-C	1 (17)	7.88 %	7.69	9
	TDC-R	1 (23)	5.70 %	9.05	6
	TDC-RC	0 (16)	10.81 %	11.31	3
	GDD-C	0 (17)	10.43 %	13.94	0
	GDD-R	6 (23)	5.36 %	8.76	13
	GDD-RC	5 (16)	8.67 %	7.36	10
	All criteria	17 (56)	4.68 %		42

Table 4 :

 4 Results for distance optimization

Acknowledgments

The authors are indebted to the anonymous referees for their numerous constructive remarks.

This research was initiated while the first author was with LAAS-CNRS, Université de Toulouse, France.