M Hinterstein

J Rouquette
email: jerome.rouquette@univ-montp2.fr

J Haines

Ph Papet

J Glaum

M Knapp

J Eckert

M Hoffman

Structural Contribution to the Ferroelectric Fatigue in Lead Zirconate Titanate (PZT) Ceramics

Many ferroelectric devices are based on doped lead zirconate titanate (PZT) ceramics with compositions near the morphotropic phase boundary (MPB), at which the relevant material's properties approach their maximum. Based on a synchrotron X-ray diffraction study of MPB PZT, bulk fatigue is unambiguosly found to arise from a less effective field induced tetragonal-tomonoclinic transformation, at which the degradation of the polarisation flipping is detected by a

less intense and more diffuse anomaly in the atomic displacement parameter of lead. The time dependance of the ferroelectric response on a structural level down to 250 μs confirms this interpretation in the time scale of the piezolectric strain response.

While the puzzle of thin film fatigue has been solved to a large degree in the last two decades , the mechanisms for bulk actuator fatigue remains still unsolved and particularly the structural contribution to the ferroelectric fatigue in bulk materials is not well understood in spite of intensive study [START_REF] Scott | Ferroelectric Memories[END_REF] .

Ferroelectric fatigue is either proposed to be caused by different mechanisms on the macroscopic or mesoscopic length scales, such as the formation of interfacial layers between the metal electrodes and the ferroelectric material [START_REF] Larsen | [END_REF]8 , 90-degree domains [START_REF] Pan | Ceramic Transactions, Ferroelectric Films[END_REF] and resultant strains, electromigration of oxygen vacancies to form extended defects capable of pinning domains 9,10 , interfacial nucleation inhibition mechanisms induced by charge injection 11,12 , macro-and/or micro-cracking 13 , etc., but to the best of our knowledge, description for fatigue on the crystal structure level has never been proposed for bulk samples. In the literature, X-ray diffraction studies of fatigue as a function of the electric field has focused on microstructural effects, which are generally probed at ψ=0° between the incident X-ray beam and the electric field vector, as the large piezoelectric strain does not permit structure refinement to be performed. PZT ceramics with compositions at the morphotropic phase boundary (MPB) [START_REF] Jaffe | Piezoelectric Ceramics[END_REF] between the tetragonal (P4mm) and rhombohedral (R3m) phase fields exhibit optimal ferroelectric properties for a large number of technological applications (sensors and actuators, micromechanical systems, and high frequency devices). We recently reported an in situ description of the origin of the ferroelectric properties as a function of the applied electric field E based on a synchrotron X-ray diffraction study 15 at ψ=45° between the incident X-ray beam and the electric field vector to minimize microstructural effects.

The monoclinic (pseudo-rhombohedral)/ tetragonal phase ratio was found to increase with electric field, which strongly supported the hypothesis of Noheda et al., who explained the strong piezoelectric properties of PZT by the presence of a monoclinic phase 16 giving rise to a polarization rotation mechanism 17 . Note that in this study as the rhombohedral model could not correctly reproduce the low-field data, this structural model was abandoned in order to describe the electric field dependence of the structural data using the same model and the controversies between the martensitic theory and the polarization rotation models are clearly beyond the scope of the present study. Additionally, polarization flipping of polar lead atoms could unambiguously be characterized by a maximum in the disorder of lead, i.e. B iso (Pb) m , which is obtained from its isotropic atomic displacement parameter linked to the maximum of the entropy for the positive-negative value of the coercive field (E c ) in the P S -E hysteresis cycle. Based on this successfull microscopic description of the macroscopic ferro-and piezolectric properties of PZT, we decided to investigate the ferroelectric fatigue in this important technological material by synchrotron X-ray diffraction. In addition, while fatigue is found to be mainly an interfacial effect in thin films, all our experiments are performed on bulk ceramics as are used in actuators, sensors and ultrasonic motors applications. Here, we unambiguosly show the structural difference between a non-fatigued and a fatigued sample which arises from i) a less effective field induced tetragonal-tomonoclinic transformation and ii) a significantly less intense and more diffuse anomaly of B iso (Pb) m , which accounts for polarisation flipping. Finally, stroboscopic pump-probe experiments are performed in order to investigate the structuralcontribution to the time dependance of the ferroelectric response down to 250 μs which clearly support and spectacularly confirm our interpretation on the structural description of the ferroelectric fatigue in the time scale of the piezolectric strain response.

X-ray measurements were performed at the high resolution powder diffraction station of the materials science beam line at the Swiss Light Source (SLS, Villigen, Switzerland) 18 using high-energy photons λ = 0.443187 Å to study thin silver metallized PZT pellets 19 (thickness ≤ 110 µm) in transmission geometry. Commercially available PIC 151 (PI Ceramics, Lederhose, Germany), which we already used [START_REF] Hinterstein | [END_REF] as it presented optimal piezo-and ferroelectric properties properties (d 33 = 500 pC/N), long operation periods and offers a high degree of reproducibility, was used to compare unfatigued and fatigued samples. The experimental set up can be found elsewhere [START_REF] Hinterstein | [END_REF] . The time to measure a complete ferroelectric hysteresis cycle was about 20 min, resulting in a cycling frequency of 0.8 mHz. The kinetics of the field induced response were investigated with a stroboscopic pump-probe setup. In this technique, the reaction is first triggered with an electric field (pump) and then a diffraction pattern is collected at a specific time delay (probe).

During a stroboscopic experiment the sample has to be cycled many times. Special care has to be taken to remain in a cycling range where no fatigue occurs. A normal measurement needs 10 s of exposure time. For bipolar cycling with 50 Hz and a time resolution of 1 ms, collection of the whole cycle yields 20 diffraction patterns. Since the setup at the MS beamline allows the measurement of only one time channel at a time, just 1/20 of the overall X-ray intensity is collected during every measurement. This results in an exposure time of 200 s per pattern and 4000 s for a whole cycle. This means that during the complete measurement, the sample runs through 4000s • 50Hz = 2•10 5 cycles. Based on the to the literature, fatigue can be detected above 10 3 cycles 20 . However, we took special care to prevent interface effects by using high quality electrodes and we clearly checked that nearly no fatigue could be observed in our material in the range of 10 5 cycles 19 . In order to assure stable conditions during the stroboscopic measurements, the samples were pre-cycled 10 5 times with a bipolar electric field of 2.0 kV/mm and a frequency of 50 Hz. In order to investigate the fatigue mechanisms, the samples were bipolar fatigued with 10 7 cycles. Fig. 1 shows the difference in the physical (Fig. 1a and1b) and structural behavior (Fig. 1c and1d) for an unfatigued and a fatigued sample. As already explained in the introduction, the ferroelectric fatigue is characterized by the loss of the switchable remanent polarization and strain as a function of the number of switching cycles. Figs. 1c and1d show the electric field dependence using a cycling frequency of 0.8 mHz of the pseudocubic 200 c reflection, which was found to be quantitatively representative of the relative abundance of the monoclinic (pseudo-rhombohedral) and tetragonal phases [START_REF] Hinterstein | [END_REF] . It is consistent with a decrease of the tetragonal intensity, while the monoclinic intensity in between the tetragonal reflections increases with electric field. However, one can clearly observe that the field induced tetragonal-to-monoclinic transformation is not as effective for the fatigued sample. This is supported by the information obtained by Rietveld refinements of the synchrotron X-ray diffraction patterns (Table 1), Fig. 2. Fig. 2a and2b shows typical high quality refinements obtained for unfatigued and fatigued PZT respectively. All the data, i.e. 24 hysteresis cycles or 2880 data sets for the unfatigued sample and 10 hysteresis cyles or 1200 data sets for the fatigued sample respectively, were refined with the same monoclinic-tetragonal two-phase model described elsewhere [START_REF] Hinterstein | [END_REF] . One can first note that the monoclinictetragonal phase fractions obtained on Fig. 3a and 3b for the unfatigued and fatigued samples respectively are consistent with a butterfly shape, which is the characteristic macroscopic signature of these materials. Additionally, the information obtained from these measurements does not exhibit dispersion and is perfectly reproducible. The edges of the butterfly shape, i.e. about 1 kV/mm, which corresponds to the increase in the monoclinic phase fraction clearly defines the coercive field as already reported [START_REF] Hinterstein | [END_REF] and this value does not change under cycling either for the unfatigued or for the fatigued samples in agreement with macroscopic measurements obtained in Fig. 1a and1b. The only large difference between the unfatigued and the fatigued samples is therefore linked to the tetragonal-tomonoclinic transformation, which is roughly shifted by 10% towards the tetragonal phase for the fatigued sample. In agreement with Noheda et al. 16 , we attributed the strong piezoelectric response of PZT to the field-induced tetragonal-to-monoclinic transition. Therefore the 10%-phase ratio shift towards the tetragonal phase for the fatigued sample is consistent with the degradation of the macroscopic properties observed on Fig. 1a and1b. Furthermore, Fig. 3c and3d show the electric field dependence of the mean square displacement of the lead atom [B iso (Pb) m ] in the monoclinic phase for the unfatigued and the fatigued sample respectively; the monoclinic phase was chosen as it is the most representative, particularly in the high electric field region. B iso (Pb) m is associated with the disorder of the considered atom 21,[START_REF] Hinterstein | [END_REF] and is therefore very sensitive to any structural change. B iso (Pb) m exhibits two anomalies at about 1kV/mm and -1kV/mm, which are consistent with the value of the coercive field. Additionally, these anomalies are always linked to a change in sign of the spontaneous polarization [( 1) and (3) in Fig. 1a] and were therefore associated to the polarization flipping of the lead atoms [START_REF] Hinterstein | [END_REF] . This means that the polarization flipping can unambiguously be detected by our structural study. One has to keep in mind that lead atoms are located in non centrosymmetric positions, i.e. they have polar displacements in the PZT structure. It is, therefore, reasonable to associate the polarization flipping, carried out by these polar atoms with a maximum in the disorder parameter of lead, i.e. Biso (Pb) that characterizes the maximum of the entropy.

Finally in order to investigate the kinetics of the field induced phase transition, we performed stroboscopic pump-probe experiments with a frequency of 50 Hz. Fig. 4 shows the characteristic 200 c reflection for the unfatigued and fatigued sample respectively as a function of the applied electric field for rectangular on-off stroboscopic pump-probe measurements. It clearly demonstrates that the field induced structural response has a time scale in the range of hundreds of microseconds as observed based on macroscopic strain measurements [22][23][24] . After 1.5 ms, no further structural changes were observed.

More importantly, it provides clear evidence for the structural contribution to the ferroelectric fatigue over a shorter timescale. The unfatigued sample exhibits a perfect electric field induced tetragonal-tomonoclinic phase transition, while the same feature is clearly incomplete for the fatigued sample.

In conclusion, in this study the ferroelectric fatigue of commercial PZT ceramics was investigated by in situ X-ray synchrotron diffraction experiments as a function of the applied electric field. Using the ψ=45° scattering geometry between the incident X-ray beam and the electric field vector previously described to minimize microstructural effects, we unambiguously characterized contribution of the crystal structure to the fatigue based on a reduced degree of tetragonal-to-monoclinic transformation, which accounts for the reduction of the piezoelectric/ferroelectric efficiency. Loss of the switchable remanent polarization could also be evidenced by a significantly less intense and more diffuse anomaly of B iso (Pb) m for the fatigued sample, which accounts for polarisation flipping. We were also able to characterize structurally the field induced response on hundred-microseconds time scale as observed based on macroscopic strain measurements 22,23 , which confirms without any doubt the structural description of the ferroelectric fatigue on a shorter timescale. These results will be useful to design and optimize high performance, fatigue free materials.

FIG. 1 .

 1 FIG. 1. Difference in the physical and structural behavior for an unfatigued and a fatigued PZT. (a)

FIG. 3 .

 3 FIG. 3. Structural dependence of PZT sample as a function of the applied electric field obtained by

Table 1 :

 1 Cell parameters, fractional atomic coordinates, isotropic (iso) atomic displacement parameters (Å 2 ) and agreement factors (%) for monoclinic and tetragonal PZT as a function of E for unfatigued sample and sample previously fatigued at 50 Hz respectively. Due to the use of hard X-ray photons, only the lead isotropic displacement parameters were refined. Cm: site occupancies Pb, Zr/ Ti, O(1) 2a (x,0,z), O(2) 4b (x, y , z). Origin at Pb position.

		Unfatigued PZT	PZT previously fatigued at 50 Hz
	Field increase		0 kV/mm		0 kV/mm
		P4mm	Cm	P4mm	Cm
	A	4.03645 (3)	5.73082 (8)	4.03641 (3)	5.72980 (12)
	B		5.72465 (10)		5.72645 (15)
	C	4.10067 (6)	4.08035 (8)	4.10052 (6)	4.08194 (10)
			90.3167 (9)		90.310 (1)
	x-Zr, Ti		0.4876 (11)		0.4852 (17)
	z-Zr, Ti	0.5414 (12)	0.5498 (11)	0.5410 (13)	0.5477 (15)
	B iso -Pb	1.335 (33)	2.37 (3)	1.15 (4)	1.82 (4)
	x-O(1)		0.4430 (38)		0.4484 (50)
	z-O(1)	0.5614 (38)	0.0323 (54)	0.5707 (34)	0.0502 (66)
	x-O(2)		0.1974 (27)		0.1880 (34)
	y-O(2)		0.2469 (23)		0.2472 (30)
	z-O(2)	0.0834 (42)	0.5864 (26)	0.0825 (45)	0.5879 (35)
	Phase Fraction	35.7 (5)	64.4 (6)	42.5 (7)	57.5 (8)
	Bragg R-factor	5.57	3.55	6.88	5.05
		Unfatigued PZT	PZT previously fatigued at 50 Hz
	Field increase	0.5694 kV/mm	0.5176 kV/mm
		P4mm	Cm	P4mm	Cm
	A	4.03666 (3)	5.73042 (8)	4.03634 (3)	5.72921 (12)
	B		5.72386 (10)		5.72493 (16)
	C	4.10154 (5)	4.08171 (10)	4.10130 (5)	4.08293 (10)
			90.3223 (8)		90.3150 (13)
	x-Zr, Ti		0.4861 (10)		0.4848 (17)
	z-Zr, Ti	0.5421 (11)	0.5486 (11)	0.5393 (13)	0.5477 (15)
	B iso -Pb	1.39 (3)	2.44 (3)	1.21 (4)	1.83 (4)
	x-O(1)		0.4429 (41)		0.4455 (53)
	z-O(1)	0.5711 (30)	0.0330 (57)	0.5722 (32)	0.0503 (68)
	x-O(2)		0.1950 (28)		0.1860 (35)
	y-O(2)		0.2436 (24)		0.2449 (32)
	z-O(2)	0.0824 (37)	0.5853 (27)	0.0783 (43)	0.5827 (36)
	Phase Fraction	39.5 (5)	60.5 (6)	44.8 (8)	55.2 (8)
	Bragg R-factor	5.43	3.87	6.76	4.89
		Unfatigued PZT	PZT previously fatigued at 50 Hz
	Field increase	0,9988 kV/mm		1 kV/mm
		P4mm	Cm	P4mm	Cm
	A	4.03728 (3)	5.73071 (10)	4.03669 (3)	5.73020 (4)
	B		5.72231 (11)		5.72154 (16)
	C	4.10164 (5)	4.08114 (8)	4.10128 (5)	4.08037 (11)
			90.3122 (10)		90.3055 (15)

This work was performed at the Swiss Light Source, Paul Scherrer Institute, Villigen, Switzerland. The authors would like to thank F. Gozzo and A. Cervellino for the beamline support and A. Bergamaschi for the competent help with the MYTHEN detector system. The research leading to these results has received funding from the BMBF (Bundesministerium für Bildung und Forschung) project grant no. 05K13PSA and from the Alexander von Humboldt Foundation project grant no. 3.1 -DEU/1150785. The authors appreciate the financial support of the German Research Foundation (DFG) through the Sonderforschungsbereich 595 'Electric fatigue in functional materials' and the German Crystallographic Society (DGK).


Mohamed Siala

Emmanuel Hebrard

Marie-José Huguet

A Study of Constraint Programming Heuristics for the Car-Sequencing Problem

Keywords: Car-sequencing, Heuristics, Tree Search Methods, Constraint Programming

), in fact, often as important or more than the propagation method. Moreover, through a criteria analysis, we are able to get several new insights into what makes a good heuristic for this problem. In particular, we show that the criterion used to select the most constrained option is critical, and the best choice is fairly reliably the "load" of an option. Similarly, branching on the type of vehicle is more efficient than branching on the use of an option. Overall, we can therefore indicate with a relatively high confidence which is the most robust strategy, or at least outline a small set of potentially best strategies.

 stating that the notion of slack used in heuristics induces a pruning rule, we propose an algorithm for this method and experimentally evaluate it, showing that, although computationally cheap and easy to implement, this is in practice a very efficient way to solve carsequencing benchmarks.

Introduction

The car-sequencing problem comes from the automotive industry and has a long history in constraint programming. In this problem, a number of cars have to be sequenced on an assembly line. Each class of cars requires a set of options. However, the working station handling a given option can only mount it on a fraction of the cars passing on the line. Each option j is then associated with a fractional number p j /q j standing for its capacity (no more than p j cars with this option j should occur in any sub-sequence of length q j ).

The car-sequencing was first introduced as a Constraint Satisfaction Problem (CSP) modeled using CHIP [START_REF] Dincbas | Solving the Car-Sequencing Problem in Constraint Logic Programming[END_REF], and later an optimization variant of this problem was used as benchmark for the Roadef Challenge in 2005 (see [START_REF] Solnon | The car sequencing problem: Overview of state-of-the-art methods and industrial case-study of the ROADEF'2005 challenge problem[END_REF] for a survey of the methods used during the Challenge). The decision version of the problem, that we shall consider in this paper, asks if a given set of cars can be sequenced on the assembly line in such a way that the line never has to slow down. This problem is NP-hard in the strong sense [START_REF] Kis | On the complexity of the car sequencing problem[END_REF].

Our contributions are divided in two major points. First, we provide an empirical study regarding car-sequencing heuristics. To the best of our knowledge, this is the most complete heuristics study for this problem. For the experiments, all heuristics are included in a complete chronological backtracking method. Moreover, we combine these heuristics to several known filtering algorithms to evaluate the trade-off between search and propagation.

We show the interest of some new heuristics in our experiments as well as the impact of the branching strategy comparatively to filtering algorithms.

This empirical study is based on a new classification of heuristics for the car sequencing. This classification is based on a set of four criteria: branching variables, exploration directions, selection of branching variables and aggregation functions for this selection. In particular, we show that the way of selecting the most constrained option is critical, and the best choice is fairly reliably the "load" of an option, that is the ratio between its demand and the capacity of the corresponding machine. Similarly, branching on the class of vehicle is more efficient than branching on the use of an option.

Overall, even though results vary greatly from instance to instance, we can therefore indicate with a relatively high confidence which is the most robust strategy, or at least outline a small set of potentially best strategies.

Second, we propose a new filtering rule, using the notion of slack introduced by Puget and Régin [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF]. This filtering can be used only when variables are explored in the lexicographic order. However, as shown by Smith [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF], this is a good strategy for this problem. Moreover, this filtering rule greatly improves the performance of our model whilst being almost free computationally.

In this paper, we are mainly interested in heuristics for solving the decision variant of the car-sequencing problem, although we shall draw a link with propagation in Section 5. However, such study provides insights on resolving optimization car-sequencing variants. In fact, there are two standard types of objective functions for the car-sequencing problem. One somehow counts the number of constraint violations. This is the objective function used for instance in the Roadef challenge 2005 [START_REF] Solnon | The car sequencing problem: Overview of state-of-the-art methods and industrial case-study of the ROADEF'2005 challenge problem[END_REF]. The other considers the delay incurred because of the violations of the capacity constraints. It was proposed by Hindi and Ploszajski [START_REF] Khalil | Formulation and solution of a selection and sequencing problem in car manufacture[END_REF] and latter on used in [START_REF] Perron | Combining Forces to Solve the Car Sequencing Problem[END_REF][START_REF] Perron | Propagation Guided Large Neighborhood Search[END_REF][START_REF] Siala | An Optimal Arc Consistency Algorithm for a Particular Case of Sequence Constraint[END_REF]. The latter objective is a simple and effective way of setting up the assembly line so that it can run smoothly. In this case, the number of extra empty slots gives an exact measure of the delay incurred. Solving a sequence of the satisfaction problem where a number of extra empty slots are added, for instance through binary search, is a natural way of optimizing this objective function. Furthermore, methods based on constraint programming are still competitive and branching heuristics are significantly contributing to this. For instance, they are relevant within a Large Neighborhood Search (LNS) approach such as the one proposed in Perron and Shaw [START_REF] Perron | Combining Forces to Solve the Car Sequencing Problem[END_REF][START_REF] Perron | Propagation Guided Large Neighborhood Search[END_REF]. In this context, one needs to solve a small CSP standing for a small part of the problem. Solving this sub-problem is expected to be fast, hence we rely on branching and selection heuristics to speed up this task. The heuristics are therefore likely to make a significant difference.

The rest of the article is organized as follows. In Section 2, we describe the car-sequencing problem and discuss the related constraint satisfaction models. In Section 3, we propose and classify a number of new and existing heuristics, which we empirically evaluate and analyze in Section 4. Then, in Section 5, we introduce a simple filtering rule to propagate capacity con-straints coupled with cardinality. Finally, we show in Section 6 the interest of the proposed filtering rule and provide a comparison with state-of-the-art propagators for this problem as well as other approaches.

The Car-sequencing problem

Problem description

In the car-sequencing problem, n vehicles have to be produced on an assembly line. There are k classes of vehicles and m types of options. Each class c ∈ {1, . . . , k} is associated with a demand D class c , that is, the number of occurrences of this class on the assembly line, and a set of options O c ⊆ {1, . . . , m}. Each option is handled by a working station able to process only a fraction of the vehicles passing on the line. The capacity of an option j is defined by two integers p j and q j , such that no subsequence of size q j may contain more than p j vehicles requiring option j.

A solution of the problem is then a sequence of cars satisfying both demand and capacity constraints.

For convenience, we shall also define, for each option j, the corresponding set of classes of vehicles requiring this option

C j = {c | j ∈ O c },
• O c 1 = {1, 3}, O c 2 = {2, 3}, O c 3 = {4}. • D class c 1 = 2, D class c 2 = 3, D class c 3 = 2
• p i /q i (lexicographically): 1/2; 2/3; 3/5; 3/6.

From above, we obtain:

• C 1 = {1}, C 2 = {2}, C 3 = {1, 2} and C 4 = {3} • D 1 = 2, D 2 = 3, D 3 = 5 and D 4 = 2 The sequence [c 2 , c 2 , c 1 , c 3 , c 3 , c 2 , c 1 ]
is a possible solution for this instance.

Constraint Programming

Many combinatorial problems can be easily formulated using a constraintbased model. In this approach, the problems are usually solved through complete search tree methods (we refer the reader to [START_REF] Marriott | Programming with constraints: an introduction[END_REF][START_REF] Rossi | Handbook of Constraint Programming[END_REF] for a comprehensive introduction). The efficiency of this approach depends on several factors: the heuristics used to select the most promising variables and values; how the search tree is expanded and how the constraints are propagated during the search.

The resolution methods that we shall study in this paper are based on depth-first search with chronological backtracking and propagation. In this context, we need to decide at each step on which variable should we branch, and which value it should be assigned to. Both decisions are based on heuristics. Such heuristics can be defined declaratively as a variable and a value ordering, respectively. We call their combination a branching strategy. We explore then the search space in a depth-first manner following the variable and value ordering prescribed by the branching strategy. The efficiency of this kind of methods clearly relies on the number of expanded nodes during search. Moreover, the shape and the size of the search tree is highly dependent on the variable and value orderings. The basic principles explaining the efficiency of variable and value orderings are well known. On the one hand, we want to choose the value that has the best chances to lead to a solution, if the current sub-tree has one. This has been called promise [START_REF] Andreas | Dual viewpoint heuristics for binary constraint satisfaction problems[END_REF] or succeed-first principle. On the other hand, we want to choose the variable that will lead to fail as soon as possible, when the current sub-tree is inconsistent.

The other crucial factor when exploring the search space is propagation. A propagator (called also filtering algorithm) associated to a constraint is a procedure that removes, for some variables, some values (or set of values) that cannot satisfy the constraint w.r.t the current assignment. The set of propagators are repeatedly called whenever a domain change occurs during search until no more reduction is possible. A good trade-off between heuristics and filtering algorithms is mandatory for better performance.

Constraint-based models for the car-sequencing

We use a standard CP modeling with two sets of variables. The first set corresponds to n integer variables {x 1 , . . . , x n } (called class variables) taking values in {1, . . . , k} and standing for the class of vehicles in each slot of the assembly line. The second set of variables corresponds to nm Boolean variables {y 1 1 , . . . , y m n } (called option variables), where y j i stands for whether the vehicle in the i th slot requires option j.

There are three sets of constraints. . This constraint is usually enforced with a Global Cardinality Constraint (Gcc) [START_REF] Charles | Generalized Arc Consistency for Global Cardinality Constraint[END_REF][START_REF] Quimper | An Efficient Bounds Consistency Algorithm for the Global Cardinality Constraint[END_REF] 2. Capacity constraints: for each option j ∈ {1..m}, no subsequence of size q j involves more than p j options of this type:

i+q j -1 l=i y j l ≤ p j , ∀i ∈ {1, . . . , nq j + 1}. In order to factor out as much as possible the propagation aspect from the study, we used several models in order to diversify the data set. More precisely, we shall consider four models, differentiated by how capacity constraints are modeled and thus propagated. For each option j, these constraints can be expressed in one of the following alternatives: (a) a naive decomposition using sum constraints (denoted by Sum). (b) a model using the Global Sequencing Constraint (Gsc) proposed by Régin and Puget [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF]. (c) a model using the AtMostSeqCard constraint recently proposed in [START_REF] Siala | An Optimal Arc Consistency Algorithm for a Particular Case of Sequence Constraint[END_REF]. (d) a model combining the Gsc and the AtMostSeqCard propagators (since the pruning obtained by both approaches is incomparable). 3. Channeling: option and class variables are channeled through simple constraints: y j i = 1 ⇔ j ∈ C x i , ∀j ∈ {1, ..., m}, ∀i ∈ {1, ..., n}. Each such constraint is implemented using a set of simple binary constraints

x i = c ⇒ y j i = 1, ∀j ∈ O c and x i = c ⇒ y j i = 0, ∀j ∈ {1, . . . , m} \ O c .

Related work

Regarding the search strategy, two main principles are known to be important for the car-sequencing problem. First, the sequence of variables to branch should follow the assembly line itself. Indeed, the structure in chain of capacity constraints makes it difficult to achieve any inference far away from a modified variable in the sequence [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF]. Second, one should assign the most constrained class or option first. This has been perceived as a fail-first strategy, hence surprising since succeed-first strategies should be better for selecting the next branch to follow. However, as pointed out in [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF], since the solutions to this problem are permutations of a multiset of values, choosing the most constrained one when it is still possible actually yields the least constrained sub-problem. In fact, in this sense, it is indeed a succeed-first strategy.

In [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF], a lexicographical exploration of the integer variables x 1 up to x n , standing for classes of vehicles, was advocated as an interesting search strategy. Three parameters were considered for choosing the most constrained class: the number of options per class (denoted as max option), the tightness of each option (ie. the capacity constraint q/p) and the usage of each option (i.e. usage rate d.q/p n ). In [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF], the authors proposed to branch on option variables y j i , exploring the sequence consistently with their position on the assembly line, however starting from the middle towards the extremities. Indeed variables at bothe ends are subject to fewer capacity constraints than variables within the sequence. Moreover, they introduced for the first time the notion of slack for selecting the most constrained option.

In [START_REF] Gottlieb | A Study of Greedy, Local Search and Ant Colony Optimization Approaches for Car Sequencing Problems[END_REF], several heuristics were compared for solving an optimization variant of this problem. These heuristics are based on the usage rate previously defined for selecting the next variables to assign in the sequence. They consider two ways for aggregating these values (using the maximum value, lexicographically, or a simple sum) when branching on class variables. Two possibilities of using the usage rate were compared : static and dynamic (i.e. updated at each node). Note that the static values of usage rate, load or slack are all equivalent. Their experiments showed essentially the interest of dynamic heuristics comparatively to static ones. The same observation is made in [START_REF] Boivin | Résolution du problème de car-sequencing à l'aide d'une approche de type FC[END_REF] where a dynamic load was used with a class variable branching and a simple summation to aggregate the values.

Heuristics Classification

Classification criteria

We propose to classify the heuristics related to this problem according to four criteria:

• The type of branching decisions: that is whether we branch on classes or options.

• The order in which we explore the variables along the assembly line: one can start from the left of the sequence and progress to the right, or start from the middle of the assembly line widen to the sides.

• The measure used to select the most constrained options.

• The function used to aggregate the evaluation of the different options in order to choose the next class of vehicles.

Notice that among the many combinations of these four criteria, some correspond to existing heuristics, however some are novel. For each criterion, there are several alternatives, we present each of them in the following.

Branching

The branching is either the assignment of a class to a slot, that is, branching on class variables x i , or the assignment of an option to a slot, that is, branching on option variables y j i . The former was used in [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF], while the latter was proposed in [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF]. Notice that when branching on option variables, we always set it to the value 1, which amounts at forcing to the corresponding option to be represented in that slot. We therefore consider these two cases denoted respectively class and opt.

Exploration

Heuristics that do not follow the sequence of variables along the assembly line generally have poor performances. Indeed, the structure in chain of capacity constraints makes it difficult to achieve any inference far away from a modified variable in the sequence [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF]. In the literature, two exploration orders were considered: either lexicographical order on class variables or from the middle to the sides of the sequence. For each type of variables, we therefore consider these two exploration cases denoted respectively lex and mid.

Selection

The best heuristics are those selecting first the most constrained option or class. Observe that since each class is defined by a set of options, then all we care about is the hardness of the options. We therefore consider the following indicators proposed in the literature to select the most constrained option:

• The capacity q j /p j : The greater the ratio q j /p j , the more constrained is the option. In fact, a greater ratio q j /p j has more impact on neighboring slots as it is shown in example 3.1.

Example 3.1. Let o 1 and o 2 be two options with 1/3 and 2/3 as capacity ratios (p j /q j ) respectively. Consider now a sequence of 5 slots in which we have to choose between o 1 and o 2 on the third position.

The two parts of the following figure show the impact of each option.

In fact, by choosing o 1 , all neighboring slots can no longer contain this option because of the at most 1/3 constraints.

y 1 i y 2 i 0 0 1 0 0 . . 1 . .
• The residual demand (d j ): This value is equal to the total demand (of a given option) minus the number of cars containing this option already allocated (d j = (D j -n i=1 min(D(y j i ))). Clearly, a greater demand makes it more difficult to fit the cars requiring this option on the line.

• The load δ j : This parameter combines the residual demand with the capacity ratio: δ j = d j q j p j . In fact, this value is tied to the number of slots required to mount d j times the option j. A greater value of the load is therefore more constrained.

• The slack σ j : Let n j be the number of slots available for option j. The slack of an option j is σ j = n jδ j . Since we want higher values to indicate more constrained options, we use in fact nσ j .

• The usage rate ρ j : This value is defined as the load divided by the number of remaining slots: ρ j = δ j /n j . It therefore represents how much of the remaining space will be occupied by vehicles requiring this option.

Based on these indicators, we consider five methods to evaluate the options. Each method returns an indicative value on how constrained is an option. In other words, the option maximizing the given parameter will be preferred in the next decision. In the following, we denote the above selection criteria respectively by q/p, d, δ, nσ and ρ. Moreover, we consider the constant function 1 as another possible selection criterion. This is proposed so that our classification also includes the max option heuristic [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF] where each class is evaluated simply by its number of options. Note also that we consider only dynamic evaluation with the four criteria : demand, load, usage rate and slack since they outperform the static versions [START_REF] Boivin | Résolution du problème de car-sequencing à l'aide d'une approche de type FC[END_REF][START_REF] Gottlieb | A Study of Greedy, Local Search and Ant Colony Optimization Approaches for Car Sequencing Problems[END_REF].

Aggregation

In the case of class branching, since classes are defined as a set of options, the decision is most often done by summing up the "scores" of the options for each class. However, there are many ways to aggregate these values. We therefore propose to add the method used for the aggregation as a fourth criterion. Let f : {1, . . . , m} → R be a scoring function. We denote f (O c ) the tuple formed by the sorted scores of class c's options, i.e., f

(O c ) = f (j 1 ), . . . , f (j |Oc| ) such that {j 1 , . . . , j |Oc| } = O c and f (j l ) ≥ f (j l+1 ) ∀l ∈ [1, . . . , |O c | -1]
. We shall consider the following ordering relations between classes:

• Sum of the elements (≤ ): c 1 ≤ c 2 iff v∈f (Oc 1 ) v ≤ v∈f (Oc 2 ) v. • Euclidean norm (≤ Euc ): c 1 ≤ Euc c 2 iff v∈f (Oc 1 ) v 2 ≤ v∈f (Oc 2 ) v 2 . • Lex order (≤ lex ): c 1 ≤ lex c 2 iff f (O c 2 ) comes lexicographically after f (O c 1 ).
Example 3.2. To illustrate aggregation functions, we consider example 2.1 and suppose that one branch on classes. In table 1, we give the different values of each selection parameter for all options. In order to emphasize the impact of aggregation functions, we propose to study the different scores for each class using the q/p parameter. Recall that each class is defined by a set of options, we obtain in table 2 the corresponding values for each class.

In table 3, we report the order of preferences given by the different aggregations. The class having the higher score will be selected first and so on.

Table 2: Classes' scores using the parameter q/p

❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ ❳ Options Classes c 1 c 2 c 3 1 2 - - 2 - 1.5 - 3 1.66 1.66 - 4 - - 2 Table 3: Scores & Heuristic decisions Agg. Scores Heuristic preferences c 1 c 2 c 3 ≤ 3.66 3.16 2 [c 1 , c 2 , c 3 ] ≤ Euc 6.75 5 4 [c 1 , c 2 , c 3 ] ≤ lex [2, 1.66, -, -] [1.66, 1.5,-,-] [2,-,-,-] [c 1 , c 3 , c 2 ]
Although we treat a simple case, one can observe that decisions can be influenced by aggregation functions. The behavior of ≤ lex here was different from the others. It prefers c 3 rather than c 2 for the second variable selection.

Heuristics structure

In the rest of the article, we denote the set of heuristics as follows: {class, opt}, {lex, mid}, {1, q/p, d, δ, nσ, ρ}, {≤ , ≤ Euc , ≤ lex } . Observe, however, that not all combinations make sense. For instance, the aggregation function does not matter when branching on options. Therefore, using the new classification, we obtain 42 possible heuristics:

• {class}, {lex, mid}, {q/p, d, δ, n-σ, ρ}, {≤ , ≤ Euc , ≤ lex } : The 30 heuristics that branches on class variables with the two exploration strategies {lex, mid}, the five selection parameters {q/p, d, δ, nσ, ρ} and the 3 aggregation techniques {≤ , ≤ Euc , ≤ lex }.

• {opt}, {lex, mid}, {q/p, d, δ, nσ, ρ}, ∅ : 10 heuristics branching on option variables with the two exploration possibilities {lex, mid} and the five selection parameters {q/p, d, δ, nσ, ρ}.

• {class}, {lex, mid}, {1}, {≤ } : The two possible heuristics related to the particular case of max option.

Among the many combinations defined by this structure, there are several existing heuristics as well as new ones. In the literature, only few heuristics have been studied.

First, the max option heuristic proposed in [START_REF] Smith | Succeed-first or Fail-first: A Case Study in Variable and Value Ordering[END_REF] branches on class variables lexicographically (lex) and the most constrained class is then selected using the sum (≤ ) aggregation. It therefore corresponds to class, lex, 1, ≤ . Second, in [START_REF] Gottlieb | A Study of Greedy, Local Search and Ant Colony Optimization Approaches for Car Sequencing Problems[END_REF], the authors proposed to use the usage rage with class branching, lexicographical exploration (lex) and ≤ , ≤ lex for aggregation. They correspond to class, lex, δ, {≤ , ≤ lex } . Similarly, the authors of [START_REF] Boivin | Résolution du problème de car-sequencing à l'aide d'une approche de type FC[END_REF] proposed a class branching using ≤ for aggregation in a lexicographical exploration (lex), however, using the load δ and the capacity q/p for selection (i.e. class, lex, {δ, q/p}, ≤ ). Finally, the heuristic proposed in [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF] is based on option branching, exploring the sequence from the middle to the sides using the slack as a selection criteria. This heuristic corresponds to opt, mid, nσ, ∅ .

To the best of our knowledge, all other heuristics are new and there is no comparative study for evaluating the impact of each classification criterion.

Evaluating the new structure

In this section, we evaluate the impact of the proposed classification criteria for the heuristics. We slightly perform randomization as a simple mechanism to deal with the Heavy tail phenomena [START_REF] Gomes | Boosting combinatorial search through randomization[END_REF]. In particular, with a low probability (2% for classes and 5% for options1 ), the second best choice (provided by the heuristic) is taken.

All the experiments were run on Intel Xeon CPUs 2.67GHz under Linux and are available via http://homepages.laas.fr/msiala/car-sequencing.

For each instance, we launched 5 randomized runs per heuristic with a 20 minutes time cut-off. All models are implemented using Ilog-Solver.

We use benchmarks available from the CSPLib [START_REF] Gent | CSPlib: a Benchmark Library for Constraints[END_REF] divided into three groups. The first group of the CSPLib contains 70 satisfiable instances having 200 cars, 5 options and from 18 to 30 classes, it is denoted by set1. The second group of the CSPLib corresponds to instances with 100 cars, 5 options and from 19 to 26 classes. In this group there are 4 satisfiable instances , denoted by set2 and 5 unsatisfiable instances denoted by set3. The third group of the CSPLib contains 30 larger instances (ranging from 200 to 400 vehicles, 5 options and from 19 to 26 classes). Set4 concerns the 7 instances from this group that are known to be satisfiable. At the top of each table, we mention, for each data set, the total number of instances with an indication on their feasibility (i.e. satisfiable: S and unsatisfiable U ). The status of the 23 remaining instances are still unknown. They are often treated in an optimization context, hence are not considered in our experimentations.

The global set of the 42 previously defined heuristics {class, opt}, {lex, mid}, {1, q/p, d, δ, nσ, ρ}, {≤ , ≤ Euc , ≤ lex } is combined with four propagators:

Sum, Gsc, AtMostSeqCard, Gsc⊕AtMostSeqCard leading to 168 different configurations. The latter is applied to each set of instances (i.e. 70 + 4 + 5 + 7 instances) with 5 randomized runs. The total average CPU time for these experiments is around 244 days.

We say that a run (related to an instance and a given configuration) is successful if either a solution was found or unsatisfiability was proven. For each set of instances, we report the percentage of successful runs (%sol )2 , the CPU time (time) in seconds both averaged over all successful runs and number of instances.

Experimental results are divided in thee parts. We first compare the many combinations of heuristic factors by giving the results for each one. Then, we study the proposed classification by evaluating each factor separately. Finally, we provide a comparison related to the efficiency and confidence of each factor

Impact of each heuristic

In this paragraph, we report the results of each heuristic separately on each set of instances averaged over the four propagators.

The set of heuristics corresponds to all possible combinations of parameters given by: {class, opt}, {lex, mid}, {1, q/p, d, δ, n-σ, ρ}, {≤ , ≤ Euc , ≤ lex } leading to the 42 heuristics presented in Section 3.

Table 4 shows the results of our experiments. For each heuristic, we indicate in column (Ref ) whether it is already known (with the corresponding reference) or not (with '-'). Recall that, in these experiments, we consider only dynamic evaluation with the four criteria : demand, load, usage rate and slack. For each set of instances, we report the percentage of successful runs (%sol) and the CPU time (time). The last two columns summarize the results over all set of instances. The column (%tot) gives the total percentage of solved instances and the column (%dev) gives the deviation in percent of a given heuristic to the heuristic solving the maximum number of instances. Bold values give the best heuristics w.r.t %sol.

For the easiest set (set1), 16 heuristics solve all instances in less than a second. Among them, 3 are known heuristics whereas 13 correspond to new combinations. It should be noted that all these configurations use a class branching and a load-based selection (i.e. ρ, δ, nσ). Interestingly, twisting one parameter from a heuristic can have a dramatic effect. For instance, the heuristic opt, lex, nσ, ∅ resolves only 32, 71% of this set whereas changing only the branching criterion to class (i.e. class, lex, nσ, {≤ lex , ≤ , ≤ Euc } ) leads to a complete resolution (i.e. 100%).

For set2 and set3, the heuristic opt, lex, ρ, ∅ gives the best results with 75% in 33.3s for set2 and 25% in 211.3s for set3. Also, the heuristics class, mid, d, {≤ lex , ≤ Euc } has the same number of successful runs compared to opt, lex, ρ, ∅ but with higher resolution time. All of these heuristics correspond to new configurations.

Finally, for set4, the best heuristics resolve 25.71% in approximately 3s and correspond to the configurations class, lex, {δ, ρ, n -σ}, ≤ lex . Another heuristic class, lex, d, ≤ lex obtains the same percentage but with higher computation time (55.3s).

Overall, the heuristic that has the best results across all data sets and therefore seems to be the more robust is class, lex, δ, ≤ lex with 86.8% of solved instances (according to the column 'Total'). More generally, heuristics using load-based selection (i.e. δ, nσ and ρ) and class branching obtain better results than the other configurations.

Criteria analysis

In this part, we aim to evaluate the relative impact of each classification criterion. For each criterion and each data set, we divide all the runs into as many sets as the number of possible values for this criterion. Then, we average the results within each set. For instance, exploration can be done either lexicographically (lex), or from the middle to the sides (mid). We will thus report two sets of statistics, one for lex and one for mid. Each average will be over one run per possible completion of the heuristic (21), filtering algorithms (4), randomized runs [START_REF] Gomes | Boosting combinatorial search through randomization[END_REF], and instances in the data set.

The following tables (5, 6, 7 and 8) are split in two parts. In the upper one, we report the results for each set and each possible criterion w.r.t the criterion being used averaged over all other criteria. The lower part shows the best results obtained for any possible combination of the other criteria. In these tables, we report the percentage of successful runs (%sol ), the CPU time (time) in seconds both averaged over all successful runs, instances and heuristic criteria. Bold values indicate best results in terms of successful runs (%sol ). Moreover, in the upper tables, the last column (%tot) gives the percentage of solved instances over the all sets.

Branching strategy

Here we compare the two branching strategies: class and opt. We tested all the possible combinations of heuristics for each strategy. However, as the constant selection parameter 1 is not defined for opt variables, we do not consider its heuristics in the following table.

When branching on opt variables, we have defined 10 heuristics (since aggregation functions are omitted): opt, {lex , mid}, {q/p, d, δ, nσ, ρ}, ∅ , that is 200 tests for each instance. To have consistent comparison with class branching, we separate its results by aggregation functions. That is class, {lex , mid}, {q/p, d, δ, n-σ, ρ}, ≤ lex , class, {lex , mid}, {q/p, d, δ, nσ, ρ}, ≤ Euc and class, {lex , mid}, {q/p, d, δ, nσ, ρ}, ≤ . The upper part of Table 5 shows that branching on classes is usually better than branching on options. However, the latter is more efficient on proving infeasibility (i.e. line opt on set3). The most efficient branching averaged over the other factors is with the aggregation ≤ but the two other aggregation (≤ lex or ≤ Euc ) are closed. This result is confirmed by the lower part of the table.

Exploration

To evaluate the exploration parameters, we consider for each ω ∈ {lex, mid} the following heuristics:

• class, ω, {q/p, d, δ, n -σ, ρ}, {≤ , ≤ Euc , ≤ lex } .
• opt, ω, {q/p, d, δ, nσ, ρ}, ∅ .

• class, ω, {1}, {≤ } . These three sets cover all possible combinations of heuristics leading to 420 tests for each parameter ω ∈ {lex, mid} and each instance. The results are shown in Table 6. In the first part of table 6, we can see that exploring the sequence from the middle then widening to the sides is in average slightly but consistently beneficial. Recall that the rationale for starting in the middle is that variables in the extremities are subject to fewer capacity constraints.

However, in the second part of table 6, we can see that in terms of successful runs, exploring the sequence using the lexicographical order leads to better results for proving unsatisfiability. This could be explained by the fact that when starting in the middle of the sequence, we effectively split the problem into essentially disjoint subproblems (there is actually a weak link through demand constraints).

Overall, the exploration parameter does not seem to be as critical as the branching parameter.

Selection

Here, we evaluate the selection criterion for choosing the most-constrained option. In this case, there are two possible sets of heuristics for each parameter ω ∈ {q/p, d, δ, nσ, ρ}:

• class, {lex, mid}, ω, {≤ , ≤ Euc , ≤ lex } • opt, {lex, mid}, ω, ∅
That is 8 heuristics for each ω combined with the 4 propagators and the 5 runs. We therefore have 160 tests for each instance (reported in table 7).

The special case of max option is presented separately at the end of Table 7 because the number of tested heuristics is different. In this case, there is only 2 heuristics class, 1, {lex, mid}, {≤ } , that is 40 tests for each instance. The upper part of Table 7 shows that using the load solves more instances in average over the all sets and for satisfiable sets (set1, set2 and set4) only. Surprisingly, the load gives better results than slack and usage rate, despite the fact that both slack and usage rate are defined using the load and the number of available slots in the variable's sequence. However the usage rate criteria seems to work better both in average and for the best results for unsatisfiable instances. Moreover, in the second part of the table, one can note that the demand obtains good results. This can be explained by the manner in which the benchmarks were generated. In fact, these instances, especially the hardest ones, are built in such way that they have a usage rate close to 1 [START_REF] Gent | CSPlib: a Benchmark Library for Constraints[END_REF]. Since the number of available slots is initially identical for all options, they also have the same (low) slack and the same (high) load. Therefore the heuristics based on these criteria (ie. load, slack and usage rate) cannot effectively discriminate values at the root of the search tree. However, recall that the load is defined as the product of the demand and the capacity. These two factors do not contribute equally, and therefore will favor different sets of options. In other words, one of them is bound to take a better decision, whilst the other is bound to take a worse one. We believe that this bias in the generation of the benchmarks explains the surprisingly good results of the demand (d) as well as the bad results of the capacity q/p along with the load, the slack and the usage rate.

Aggregation

Aggregation functions are only used with class branching. For each parameter ω ∈ {≤ lex ≤ ≤ Euc }, we have the 10 following heuristics combined with the propagators and the random runs (i.e. 200 tests for each ω and each instance):

• class, {lex, mid}, {q/p, d, δ, nσ, ρ}, ω

The constant parameter for selection 1 is not considered in these experiments since it is only defined with the ≤ aggregation. The results are given in Table 8. As we can see in the first part of this table, the three aggregation functions provide in average similar results except for the hardest instances (set4) where ≤ lex solved more instances. Considering all instances, ≤ solves the large number of problems. No solution were found for unsatisfiable instances as in our case, only opt branching can solve these instances (i.e. which by default do not use any aggregation function). However, regarding the best results in the second part of the table, when using ≤ lex and ≤ Euc , one can obtain better performances in terms of resolved instances.

Criteria Analysis Conclusions

We have previously evaluated the average best choice of each criterion (in terms of solved instances). However, this choice is not the best on each set of instances. Instead, we can determine the best choice for each data set, called the "perfect" choice. The Confidence of the average best choice can then be defined by the ratio between the average best choice and the perfect choice. Similarly, we can consider the "worst" choice for each data set, and subsequently, define the Significance of a given factor using the ratio between the worst and the perfect choice as 1worst/perf ect. In Table 9, we give the values of Confidence and Significance for each factor (branching, selection, exploration and aggregation)This table shows that there is high confidence for each selected average best choice (between 0.989 and 1.0): that is, exploration from middle to sides using a class branching, load selection, and a sum aggregation. When considering the significance of each criterion, one can observe that only two of them (branching and selection) have a valuable impact. For the two other criteria (i.e. exploration and aggregation), there is little impact on the results when changing the parameters.

Therefore, the most robust heuristics will be those branching on classes variables and selecting options using the load criterion, that is class, {lex, mid}, δ, {≤ , ≤ Euc , ≤ lex } .

Slack-based Filtering Rule

When analyzing the heuristics, we have seen that selecting the options using the load, the slack, or the usage rate is beneficial. In this section, we shall see that one can go one step further, and use the same idea to prune the search tree at a very cheap computational cost.

Slack and Pruning the Search Tree.

In [START_REF] Régin | A Filtering Algorithm for Global Sequencing Constraints[END_REF], it is observed that if the slack (σ j ) of an option j is negative, then the problem is unsatisfiable. Indeed, the load (δ j ) tends to represent the number of required slots to mount all the occurrences of an option. Since the slack is the difference between the available number of slots and the load, a negative value suggests infeasibility since we need more slots than are available. However, one has to be careful about boundaries issues since the capacity constraints are truncated at the extremities of the assembly line. For instance, consider an option j with p j = 1, q j = 3 and d j = 2. The slack is negative as soon as there are less then six slots remaining (n j < 6), however a line with only four slots is sufficient if we put the two classes requiring this option on both ends of the line. In other words, the load is an accurate measure of how many slots are needed for a given option, however only for large values of demand and length of the assembly line.

In order to take into account this issue of boundaries, we propose the following alternative definition for the load.

δ ′ j = q j (⌈d j /p j ⌉ -1) + p j if d j mod p j = 0 d j mod p j otherwise Notice however that the formula above is true only if the unassigned slots are contiguous in the assembly line. Therefore, in the following we assume that we explore the assembly line from left to right, and we describe a simple rule to filter the domain of a sequence of Boolean variables y 1 , . . . , y n subject to capacity constraints with a fixed cardinality. (i.e. q l=1 y i+l ≤ p ∀i ∈ [1, . . . , nq + 1] and n i=1 y i = d) Proposition 5.1. For each option j, δ ′ j is a lower bound on the number of required slots.

Proof. Consider a sequence of δ ′ j unassigned Boolean variables subject to a capacity constraint p j /q j and a demand d j . Now assign the first p j variables If δ ′ = ni + 1, then if d mod p = 0, we impose y i = 1 for all i such that i mod q < p. Otherwise (i.e. d mod p = 0), we impose y i = 1 for all i such that i mod q < (d mod p).

Proof. Suppose that (d mod p = 0). Then there exists two integers k and r such that d = k.p +r. Notice that in this case, we have δ ′ = q.k +r. Consider a subsequence y a , . . . , y b such that a mod q = 0 and b = a + r -1, i.e., such that the rule above applies. Then there exist two integers α and β such that

a = α • q and n -i -b = β • q (since n -i + 1 = δ ′ = q.k + r). Now using n -i -b = β • q, we show that n -i + 1 = β • q + a + r then n -i + 1 = (α + β) • q + r and hence k = α + β (since n -i + 1 = q.k + r).
However, by definition of α and β, we may argue that the number of occurrences of the value 1 on y 0 , . . . , y a-1 is at most α • p and at most β • p on y b+1 , . . . , y n-i . Now since the demand d = (α + β).p + r then all the p variables the subsequence y a , . . . , y b must take the value 1.

We use a similar argument for the second case. Suppose that d mod p = 0, consider a subsequence y a , . . . , y b such that a mod q = 0 and b = a + p -1.

Then there exist two integers α and β such that a = α•q and n-i-b = β•q. Therefore, the number of occurrences of the value 1 on y 0 , . . . , y a-1 is at most α • p and at most β • p on y b+1 , . . . , y n-i . Now using the demand d = k • p, and δ ′ = q (⌈d/p ⌉ -1) + p we show that ni + 1 = q(k -1) + p. However, since b = a + p -1, a = α • q and nib = β • q, then k = α + β + 1 and all p variables the subsequence y a , . . . , y b must take the value 1.

Figure 2 and 3 depict the proposed pruning. On the one hand, when d mod p = 0, the only possible arrangement of vehicles that satisfy the capacity constraint is to start the sequence with p vehicles requiring the option, then qp vehicles not requiring the option and repeat (see Figure 2). Notice that because of the capacity constraint, all other variables must take the value 0. On the other hand, when d mod p = 0, one must start the sequence with d mod p vehicles requiring the option, then the following q -(d mod p) slots can be filled arbitrarily as long as exactly p vehicles requiring this options are fitted in the q first slots. Here again, the initial sequence must be repeated throughout (see Figure 3).

Time Complexity

Observe that this rule is extremely cheap to enforce. Once one has computed the load, the domain filtering can be achieved in O(k) where k is the number of option variables forced to take the value 1. Indeed, when d mod p = 0 we can jump over the variables which are not forced to take the value 1, since their position is given by a simple recursion. In the worst case, i.e., when d mod p = 0, this is the same time complexity as achieving arc consistency on the AtMostSeqCard constraint [START_REF] Siala | An Optimal Arc Consistency Algorithm for a Particular Case of Sequence Constraint[END_REF]. However whilst the propagator of AtMostSeqCard always achieves its worst case complexity, this rule rarely does in practice.

Evaluating the filtering rules

In order to evaluate the slack-based filtering rule, we propose to compare its results against the other propagators. Notice first that, as we mentioned in section 5, this rule can be applied only with lex branching. In fact, we can use the following set of heuristics {class, opt}, lex, {1, q/p, d, δ, nσ, ρ}, {≤ , ≤ Euc , ≤ lex } . That is 21 different heuristics for each filtering algorithm. The experiments concern 9030 configuration per propagator. Table 10 shows that the extra filtering of the rule we introduced or from the existing ones (i.e. AtMostSeqCard, Gsc and Gsc⊕ AtMostSeqCard) does help a lot. For instance, at least 90% of the instances of the first set are resolved irrespectively of the heuristic being used against 75,89% with the default decomposition (i.e. sum). The difference is even greater for the other sets. At the outset, the gain of the new propagator seems comparable with the others. However, this method is in fact very different from the Gsc and quite close to the AtMostSeqCard constraint.

The Gsc constraint saves many more backtracks than the others. It does not subsume the pruning of our filtering rule, but it is much stronger in terms of search tree size. However the overhead of our method is negligible, whereas Gcs greatly slows down the search at the same time that it reduces its size.

Consider now the propagation method as a fifth criterion (i.e. in addition to the heuristic factors). We calculated its Confidence and Significance according to the same formula given in Section 4.3. Their values are equal to (respectively) 0.996 and 0.217. This is similar to all other criteria in terms of confidence (i.e. between 0.989 and 1.0), but slightly less than the Significance of branching and selection. This emphasizes the importance of these factors which are at least as important as the propagation level.

Overall, we observe that the choice of the search strategy has a very significant impact on the efficiency of the method. For instance, on the set of easiest instances, when averaging across all heuristics, the "worst" filtering method (decomposition into sum constraints) is successful in about 20% less runs than the best (Gsc+ AtMostSeqCard).

However, now averaging across all four models, the worst heuristic opt, lex, nσ, -, is successful 56% less runs than one of the many heuristics solving all easy instances (see table 4). For harder instances (set2, 3 and 4), these choices are even more important, with a 42% gap between the best and worst model, whilst the worst heuristics (in this case opt, lex, p/q, -) do not solve any instances. It is hardly a surprise to observe that the choice of search strategy is a critical one. However, whilst the aim of this study was to better understand what makes a good heuristic for the car-sequencing problem, it was relatively surprising to find out that minor variations around known heuristics would bring such a substantial gain.

Finally, as a comparison with other CP-based methods [START_REF] Van Hoeve | New Filtering Algorithms for Combinations of Among Constraints[END_REF][START_REF] Perron | Combining Forces to Solve the Car Sequencing Problem[END_REF][START_REF] Gottlieb | A Study of Greedy, Local Search and Ant Colony Optimization Approaches for Car Sequencing Problems[END_REF]. The first set was totally resolved in less than a second unlike other methods. The second set results are as good as the best ones. However, the Regular constraint [START_REF] Van Hoeve | New Filtering Algorithms for Combinations of Among Constraints[END_REF] maintains best results for proving infeasibility by a difference of one instance. We are not aware of any approach treating the fourth set of as decision problems. These instances are often treated in an optimization context.

Conclusion

Throughout this paper, we empirically studied a large set of heuristics for the car-sequencing problem and proposed to classify these heuristics using 4 criteria: branching variables, exploration directions, parameters for the selection of branching variables, and aggregation functions for these criteria. The experiments show the interest of some classification criteria (the branching and the selection) and a low impact of the other criteria (exploration directions and aggregation functions). Moreover, the study shows that one criterion can drastically effect the heuristic behavior.

The second contribution of this paper is the slack-based filtering algorithm. Even though quite simple and easy to implement, the filtering algorithm introduced in this paper is often as good as state-of-the-art propagators for the car-sequencing problem, and better in some cases.

A natural extension of this analysis is to study the impact of these heuristics for solving optimization variants for the car-sequencing problem. Such study gives insights for the industrial and optimization communities working on this problem. This can be achieved by decomposing the problem into a sequence of satisfaction problems or by using branching strategies directly within meta-heuristics such as Local Search and Ant Colony Optimization (see for instance [START_REF] Gottlieb | A Study of Greedy, Local Search and Ant Colony Optimization Approaches for Car Sequencing Problems[END_REF]).

.Example 2 . 1 .

 21 and the option's demand D j = c∈C j D class c Consider the simple case of 7 slots (i.e. n = 7) with 3 classes {c 1 , c 2 , c 3 } and 4 options such that:

1 .

 1 Demand constraints: for each class c ∈ {1..k}, |{i | x i = c}| = D class c

Figure 2 :

 2 Figure 2: Filtering when d mod p = 0

Figure 3 :

 3 Figure 3: Filtering when d mod p = 0

Table 1 :

 1 Values of the selection criteria for each option

	❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ Selection parameter Options	1	2	3	4
	1	1	1	1	1
	q/p	2	1.5 1.66	2
	d	2	3	5	2
	δ	4	4.5 8.33	4
	n -σ	4	4.5 8.33	4
	ρ	0.57 0.64 1.19 0.57

Table 4 :

 4 Comparison of heuristics averaged over propagation rules

		Heuristics		Ref.			Instances		Total
						set 1 (70, S) set 2 (4, S)	set 3 (5, U )		set 4 (7, S)
	Sel. Br. Expl. Aggr.		%sol time %sol time %sol time		%sol time %tot %dev
				≤ lex	[6]	100.00 0.6 52.50 59.1	0.00	-25.71	2.9 85.93 1.00
			lex	≤	[6]	100.00 0.6 48.75	0.2	0.00	-10.71 84.4 84.53 2.61
	ρ	class	mid	≤ Euc ≤ lex ≤	---	100.00 0.6 30.00 99.92 0.5 53.75 163.5 0.2 100.00 0.5 51.25 236.6	0.00 0.00 0.00	-12.85 156.3 83.84 3.42 -16.42 50.0 85.17 1.88 -18.57 5.4 85.29 1.74
				≤ Euc	-	100.00 0.5 51.25 249.3	0.00	-17.14 30.2 85.17 1.88
		opt	lex mid	--	--	87.00 1.9 75.00 33.3 25.00 211.3 87.64 2.9 31.25 0.4 23.00 233.6 14.28 171.1 75.29 13.26 5.71 533.4 76.22 12.19
				≤ lex	-	100.00 0.6 52.50 59.2	0.00	-25.71 2.8 85.93 1.00
			lex	≤	-	100.00 0.6 48.75	0.2	0.00	-10.71 78.6 84.53 2.61
	n -σ	class	mid	≤ Euc ≤ lex ≤	---	100.00 0.6 48.75 100.00 0.6 53.75 169.7 0.1 100.00 0.5 51.25 236.9	0.00 0.00 0.00	-10.71 79.4 84.53 2.61 -18.57 33.1 85.41 1.61 -22.14 29.0 85.58 1.41
				≤ Euc	-	99.92 0.5 51.25 236.3	0.00	-22.14 28.8 85.52 1.48
		opt	lex mid	--	-[14]	32.71 21.7 43.75 236.8 13.00 190.7 38.14 13.0 26.25 33.7 18.00 260.8	0.00 0.00	-29.42 66.11 -33.31 61.62
				≤ lex	-	100.00 0.6 71.25 42.4	0.00	-25.71	3.0 86.80 0.00
			lex	≤	[1]	100.00 0.6 48.75	0.3	0.00	-10.71 100.2 84.53 2.61
	δ	class	mid	≤ Euc ≤ lex ≤	---	100.00 0.6 48.75 100.00 0.5 37.50 38.2 0.3 100.00 0.5 68.75 167.9	0.00 0.00 0.00	-10.71 87.3 84.53 2.61 -15.00 51.5 84.36 2.81 -20.71 42.8 86.28 0.60
				≤ Euc	-	100.00 0.5 68.75 166.5	0.00	-20.00 16.2 86.22 0.67
		opt	lex mid	--	--	98.57 1.2 36.25 111.7 98.92 3.7 43.75 3.8	0.00 0.00	-22.85 -21.42 88.8 84.29 2.89 5.8 83.78 3.48
				≤ lex	-	82.85 7.8	0.00	-	0.00	-	0.00	-67.44 22.31
			lex	≤	[1]	83.35 10.1 18.75	0.1	0.00	-	0.00	-68.72 20.84
	q/p	class	mid	≤ Euc ≤ lex ≤	---	83.42 11.3 18.75 0.09 84.71 7.9 18.75 95.7 85.35 7.7 18.75 100.9	0.00 0.00 0.00	---	0.00 0.00 0.00	-68.77 20.77 -69.82 19.56 -70.34 18.96
				≤ Euc	-	84.64 7.5 18.75 96.0	0.00	-	0.00	-69.77 19.63
		opt	lex mid	--	--	65.71 73.3 70.71 29.8 12.50 606.4 0.00 -	0.00 0.00	--	0.00 0.00	-53.48 38.38 -58.14 33.02
				≤ lex	-	90.92 1.2 37.50 47.4	0.00	-25.71 55.3 77.84 10.32
			lex	≤	-	95.07 1.9 41.25 48.5	0.00	-17.14 21.5 80.70 7.03
	d	class	mid	≤ Euc ≤ lex ≤	---	94.50 0.7 43.75 106.5 90.64 1.9 75.00 83.4 94.71 0.6 67.50 68.9	0.00 0.00 0.00	-23.57 40.2 80.87 6.83 -24.28 5.3 79.24 8.71 -13.57 53.9 81.33 6.30
				≤ Euc	-	94.57 0.6 75.00 83.2	0.00	-15.71 50.7 81.74 5.83
		opt	lex mid	--	--	73.78 2.9 56.25 79.5 77.28 13.7 43.75 5.2	0.00 0.00	--	0.71 282.0 62.73 27.73 7.85 16.5 65.58 24.45
	1	class	lex ≤ mid ≤	[17, 1] -	86.92 13.2 18.75 89.92 8.3 63.75 20.3 0.1	0.00 0.00	--	0.00 0.00	-71.62 17.49 -76.16 12.26

Table 5 :

 5 Evaluation of the branching variants

	Av. Bran.	set1 (70, S)		set2 (4, S)		set3 (5, U )	set4 (7, S)	Global
	(×200)	%sol avg bts time	%sol avg bts time %sol avg bts time %sol avg bts time %tot
	opt	73.0 102023.9 14.1	36.8 287139.5 82.0 7.9 53275.4 225.6	7.2 207502.8 107.9 62.2
	class, ≤ lex	94.9 26120.0 2.0	45.2 481410.8 84.9	0.0	-	-17.7 98707.8 22.5 80.7
	class, ≤	95.8 27209.1 2.1 46.3 327601.5 95.7	0.0	-	-12.4 156300.3 44.6 81.1
	class, ≤ Euc	95.7 27563.3 2.1	45.5 463196.6 107.9	0.0	-	-13.2 107599.7 52.9 81.0
	Best Bran.
	opt	100.0 98577.4 10.3	75.0	7251.3	0.5 40.0 46211.8 261.8 25.7 629016.8 130.7
	class, ≤ lex 100.0	184.7 0.0 100.0 730687.4 89.5	0.0	-	-28.5 29632.6 58.5
	class, ≤	100.0	184.2 0.0	95.0 904739.2 96.3	0.0	-	-25.7 34705.3 54.8
	class, ≤ Euc 100.0	184.4 0.0 100.0 211830.5 128.8	0.0	-	-28.5 47435.1 75.4

Table 6 :

 6 Evaluation of the exploration variants

	Av. Expl.	set1 (70, S)	set2 (4, S)	set3 (5, U )	set4 (7, S)	Global
	(×420)	%sol avg bts time %sol avg bts time %sol avg bts time %sol avg bts time %tot
	lex	89.2 50617.6 5.6	40.0 259229.0 46.6	1.8 52295.1 204.3 11.3 120652.6 54.2 75.5
	mid	90.3 42167.0 4.1 46.7 479360.9 126.5 1.9 54184.0 245.5 12.7 139829.4 42.8 76.8
	Best Expl.
	lex	100.0	184.8 0.0 100.0 730687.4 89.5 40.0 46211.8 261.9 28.5 29632.6 58.5
	mid	100.0	183.5 0.0 100.0 213028.8 129.1 36.0 63984.8 307.6 28.5	1357.4 9.2

Table 7 :

 7 Evaluation of the selection variants

	Av. Selec.	set1 (70, S)		set2 (4, S)		set3 (5, U )			set4 (7, S)		Global
	(×160)	%sol avg bts time %sol avg bts time %sol avg bts time %sol avg bts time %tot
	ρ	96.8	1628.8 1.0	49.2 480035.3 99.9 6.0 49922.8 222.0 15.1 136850.7 81.7	82.6
	n -σ	83.8	5773.4 2.3	47.0 699885.1 126.9	3.8 58466.5 231.4 13.7 103897.2 33.3	71.7
	δ	99.6	3292.6 1.0	52.9 254264.1 74.8	0.0	-	-18.3 98161.0 41.5	85.1
	q/p	80.0 195896.5 17.7	13.2 135511.2 123.0	0.0	-	-	0.0	-	-	65.8
	d	88.9 25988.2 2.7 55.0 254347.0 68.8	0.0	-	-16.0 185381.6 36.8	76.2
	1 (×40)	88.4 130722.2 10.7	41.2	28165.2 15.8	0.0	-	-	0.0	-	-	73.8
	Best Selec.
	ρ	100.0	184.8 0.0	75.0	7251.3	0.5 40.0 46211.8 261.9 25.7	4843.0 0.4
	n -σ	100.0	184.8 0.0	75.0 1009607.4 124.1 32.0 75445.9 351.0 25.7	4843.0 0.4
	δ	100.0	184.8 0.1 100.0 730687.4 89.5	0.0	-	-25.7	4843.0 0.4
	q/p	98.8	7208.4 3.4	25.0	68.2	0.1	0.0	-	-	0.0	-	-
	d	100.0	178.7 1.2 100.0 213028.8 129.1	0.0	-	-28.5 29632.6 58.5
	1	99.7 58773.0 9.9	85.0	51740.9 36.9	0.0	-	-	0.0	-	-

Table 8 :

 8 Evaluation of the aggregation variants

	Av. Agg.	set1 (70, S)	set2 (4, S)	set3 (5, U )	set4 (7, S)	Global
	(×200)	%sol avg bts time	%sol avg bts time %sol avg bts time %sol avg bts time %tot
	≤ lex	94.9 26120.0 2.0	45.2 481410.8 84.9	0.0	-	-17.7 98707.8 22.5	80.7
	≤	95.8 27209.1 2.1 46.3 327601.5 95.7	0.0	-	-12.4 156300.3 44.6	81.1
	≤ Euc	95.7 27563.3 2.1	45.5 463196.6 107.9	0.0	-	-13.2 107599.7 52.9	81.0
	Best Agg.
	≤ lex	100.0	184.7 0.0 100.0 730687.4 89.5	0.0	-	-28.5 29632.6 58.5
	≤	100.0	184.2 0.0	95.0 904739.2 96.3	0.0	-	-25.7 34705.3 54.8
	≤ Euc	100.0	184.4 0.0 100.0 211830.5 128.8	0.0	-	-28.5 47435.1 75.4

Table 9 :

 9 Confidence and significance for each factor

		Confidence Significance
	Branching	0.989	0.247
	Selection	0.995	0.231
	Exploration	1.000	0.017
	Aggregation	0.995	0.015

Table 10 :

 10 Evaluation of the filtering variants (averaged over all heuristics)

	Filtering (×21)	set1 (70 × 5) %sol avg bts time %sol avg bts time %sol avg bts time %sol avg bts time %tot set2 (4 × 5) set3 (5 × 5) set4 (7 × 5) Global
	sum	75.8 190636.0 11.2 22.6 792179.8 44.4	0.0	-	-	7.7 194651.7 17.0	63.4
	Gsc	94.8	1639.4 4.2 44.0	38673.7 49.2 2.8 49417.9 260.8 12.1	35302.0 64.3	80.4
	AtMostSeqCard	91.2 36285.7 3.9 49.2 411514.8 46.2	1.5 68873.9 15.1 13.1 239317.8 41.4	77.7
	Gsc ⊕ AtMostSeqCard 95.1 1585.1 4.3 44.0	35711.3 45.4 2.8 46330.2 248.6 12.5	32258.4 80.9	80.6
	Slack-based	90.5 55384.8 3.8 43.3 627443.4 43.9	1.7 82815.9 16.1 12.2 356073.4 34.8	76.7

Table 11 :

 11 Best results for filtering variants

	Filtering	set1 (70 × 5) %sol avg bts time %sol	set2 (4 × 5) avg bts time %sol avg bts time %sol avg bts time set3 (5 × 5) set4 (7 × 5)
	sum	100.0 184.8 0.0	75.0	7251.3	0.5	0.0	-	-25.7	4843.0	0.4
	Gsc	100.0	184.8	1.2	75.0	18073.7 58.2 40.0 46211.8 261.9 28.5 29632.6 58.5
	AtMostSeqCard	100.0 184.8 0.0 100.0 730687.4 89.5 20.0 60460.4	13.5 28.5 31617.6 6.0
	Gsc⊕AtMostSeqCard 100.0	184.8	1.2	75.0	16923.7 55.0 40.0 46196.7 259.7 28.5 17252.6 40.8
	Slack-based	100.0 184.3 0.0	75.0 510189.0 35.1 20.0 70573.6	14.0 28.5 332430.9 34.3

Those values were arbitrarily chosen. The impact of branching on an option variable being lower, a higher probability was necessary.

Since set3 contains only unsatisfiable instances, then %sol corresponds to the percentage of instances for which the heuristics prove the unsatisfiability

Acknowledgement

We would like to thank our Cork Constraint Computation Center colleagues for giving us access to their computing resources.

y j 0 1 1 1 2 1 3 0 4 0 5 1 6 1 7 1 8 0 9 0 10 1 11 1 12 1 13 0 14 0 15 1 16 1 3 2 6

to 1, then the q jp j next variables to 0 and repeat this (⌈d j /p j ⌉ -1) times.

Then fill the remaining variables with the value 1. The sequence built in this way is of length δ ′ j and cardinality d j . Moreover, every subsequence of length q j has exactly p j times the value 1, therefore, it is not possible to obtain the same cardinality in a shorter sequence, hence δ ′ j is a lower bound.

Filtering the Domains

We suppose now that all variables up to a rank i -1 are ground. To make the notation lighter we rename the sequence of unassigned variables y i , . . . , y n to: y 0 , . . . , y n-i .

When the real load δ ′ is greater than the residual number of slots n-i+1, then we should fail since δ ′ is a lower bound on the number of required slots. When δ ′ = ni + 1, we can filter out some values. Moreover, we can prune inconsistent values in the domains of the option variables when the load is equal to the remaining number of slots. We illustrate this situation in Example 5.1

Example 5.1. Figure 1 shows a possible assignment for a sequence y j 0 , . . . , y j 16 , with capacity 3/5 and demand 11. Consider the two slots indexed 5 and 6, corresponding to the variables y j 5 and y j 6 . On the left, there are 5 slots, hence we can fit at most 3 vehicles with the option j, indeed fitting 4 vehicles requires 6 = 5(⌈4/3⌉ -1) + 4 mod 3 slots. Similarly, on the right, one cannot fit more than 6 vehicles with option j since fitting 7 vehicles would require 11 slots. Therefore, since the total demand is 11, we can conclude that 11 -6 -3 = 2 vehicles with option j must fit in the slots 5 and 6. In other words, both y j 5 and y j 6 must be equal to 1. Now we formally define a pruning rule that can detect all such forced assignments (e.g., it detects all bold faced 1's in Figure 1). Theorem 5.1. The following filtering rule is correct:


Benjamin Ribba
email: benjamin.ribba@inria.fr

France Mentré

On The Use of Model-Based Tumor-Size Metrics To Predict Survival

Mixed

)n all these studies, the authors used a mixed-effects model to analyze the observed tumor size data. They then estimated individual patient parameters using all available tumor size measurements so called individual empirical Bayes estimates . These individual parameter estimates were used to calculate certain tumor size response metrics, which were subsequently tested as predictors of overall survival.

  -effect models are increasingly being used to analyze the time-course of tumor size and to identify tumor size metrics as predictors of overall survival in cancer patients. However, caution is needed interpreting the results of such analyses and applying them to clinical research. Sparse individual tumor data as a result of patient mortality can lead to falsely conclusions about tumor size metrics as predictors of overall survival. Body text Using early tumor size response to predict long-term clinical outcome, such as overall survival, may shorten the time needed to evaluate the effectiveness of treatments [ ]. A recent review has summarized efforts to model continuous longitudinal tumor size data to describe the observed time course of tumor size and the extension to prediction of survival [ ]. By using time-to-event approaches, metrics derived from the time course of individual tumor size changes have been tested as predictors of overall survival. Clinical Pharmacology and Therapeutics was the first journal to publish such an analysis [ ]. Similar studies have also been reported for colorectal cancer, thyroid cancer, breast cancer and gastrointestinal stromal tumor see [ ] .

  Figure .


Dima Abi Abdallah

Vincent Robin

Agnès Drochon

I. INTRODUCTION

AGNETIC Resonance Imaging (MRI) has become by far the primary tool for gaining important insights into the functional and metabolic bases of heart disease. However, MRI observations of a moving organ, such as the heart, require synchronization: since an image cannot be wholly acquired in one heart cycle, its successive acquisitions have to be accurately combined with the cardiac phase motion in order to overcome movement related blurring effects. Cardiac gating is usually done by detecting the R peaks on the simultaneously recorded electrocardiogram (ECG), which are then used to trigger the consecutive image acquisition sequences. Hence, achieving good image quality depends greatly on a reliable trigger signal in order to guarantee that accurate data collection always starts at the same point of the cardiac cycle. However, accomplishing an accurate R peak detection, and therefore a correct synchronization, is often obstructed by the high corruption levels of the ECG signal due to electromagnetic effects, especially when high resolution imaging is performed where high static magnetic fields as D. Abi-Abdallah is with the Laboratory of Biomechanics and Biomedical Engineering of the University of Technology of Compiègne, France (corresponding author e-mail: dima.abiabdallah@utc.fr).

V. Robin is with the Laboratory of Applied Mathematics of the University of Technology of Compiègne, France (e-mail: vincent.robin@utc.fr).

A. Drochon is with the Laboratory of Biomechanics and Biomedical Engineering of the University of Technology of Compiègne, France (email: agnes.drochon@utc.fr).

O. Fokapu is with the Laboratory of Biomechanics and Biomedical Engineering of the University of Technology of Compiègne, and the University of Picardie Jules Verne, France (e-mail: odette.fokapu@utc.fr).

. well as strong and fast switching gradients are needed. Most of the previous studies [START_REF] Abächerli | Suppression of MR gradient artefacts on electrophysiological signals based on an adaptive real-time filter with LMS coefficient updates[END_REF], [START_REF] Damji | Interference suppression in cardiac synchronization system operating in a high magnetic field NMR imaging system[END_REF], [START_REF] Felblinger | Restoration of Electrophysiological Signals Distorted by Inductive Effects of Magnetic Field Gradients during MR Sequences[END_REF] that developed methods to produce a synchronization signal aimed mainly to reduce gradient switching, radiofrequency pulses and patient movements generated artifacts and have not been proven to be efficient for the magnetohydrodynamic (MHD) artifact removal. To our knowledge the only solution for eliminating this artifact was proposed by Fischer et.al. [START_REF] Fischer | Novel real-time Rwave detection algorithm based on the vectorcardiogram for accurate gated magnetic resonance acquisitions[END_REF]. Their method exploits a vectocardiogram (VCG) signal and necessitates the use of multiple sensors. Nevertheless, positioning a large number of sensors in a NMR environment is not always simple, and might have consequences on the image quality. The MHD effect is of physiological origins and is due to the blood motion in the magnetic environment. In fact, the blood's charged particles flowing in the static magnetic field get deflected by the Lorentz force and as a result a Hall potential is generated across the vessel walls. The magnitude of these voltages is a function of the blood flow rate and vessel diameter, their amplitude is maximal when the blood flow is perpendicular to the magnetic field. These superimposed voltages, induced by the electrodynamic interactions of the static field with arterial blood flow, especially in the ascending aorta and aortic arch, alter the waveform of the recorded ECG. The largest magnetically induced voltages occur during the blood ejection phase in the aorta which coincides with the ventricle repolarization (T wave). Hence the major MHD-caused change observed on the simultaneously recorded ECG is an increase of the T wave amplitude during in vivo MRI observations. Magnetic field interactions with blow flow have been demonstrated by multiple authors throughout in vitro experiments [START_REF] Keltner | Magnetohydrodynamics of blood flow[END_REF] as well as in vivo studies [START_REF] Gaffey | Alterations in the rat electrocardiogram induced by stationary magnetic fields[END_REF], [START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF] where ECG alterations have been observed. For small animals, such as rats, the T wave elevation could be quantitatively measured at field levels above 0.3 Tesla , at a field intensity of 1.5 Tesla Gaffey et.al. [START_REF] Gaffey | Alterations in the rat electrocardiogram induced by stationary magnetic fields[END_REF] measured a 216% T wave elevation. However for larger species such as monkeys, who have larger aortas and greater flow rates, the MHD effect could be noted for much lower field intensities starting at 0.1 Tesla, Tenforde et.al. [START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF] measured elevations greater than 300% on macaca monkeys exposed to a 1.5 Tesla magnet. Nowadays, as field intensities augment, in the strive for achieving better image quality, the MHD effect is expected to amplify. Chakeres et.al. [START_REF] Chakeres | Effect of stati magnetic field exposure of up to 8 Tesla on sequential human vital sign measurements[END_REF] recently conducted a study with an 8 Tesla magnet to examine the influence of high magnetic fields on human vital signs; they

Alterations in human ECG due to the MagnetoHydroDynamic effect: A method for accurate R peak detection in the presence of high MHD artifacts

Dima Abi-Abdallah, Vincent Robin, Agnès Drochon, and Odette Fokapu M noted the same effects on the ECG as Gaffey and Tenforde but did not make any T wave elevation measurements.

In the cardiac MRI gating context, the MHD effect can be particularly hindering for synchronization, as the T wave might reach amplitudes as high as the QRS complex and could be mistaken for an R peak, hence leading to incorrect acquisition triggering.

In this study we examine the MHD effect on several ECG signals recorded on different subjects while exposed to a 1.5 Tesla static magnet and assess the major alterations imposed on the signal. We also review a previously elaborated real time R peak detection algorithm [START_REF] Abi-Abdallah | Cardiac and respiratory MRI gating using combined wavelet sub-band decomposition and adaptive filtering[END_REF] and test its capacity of extracting correct triggers from these MHD contaminated signals.

II. METHODS

A. Signal acquisition

The ECG signals were recorded on 6 healthy subjects aged between 20 and 40 years. MR compatible electrodes (3M Red Dot) were placed on the subjects' chest area according to the Einthoven triangle configuration, and three ECG leads signals were acquired using an optical nonmagnetic amplifier (0.5Hz-20Hz). For each subject an ECG signal was first recorded away from the magnet (at 0 Tesla) afterwards the subject was placed inside a 1.5 Tesla Signa GE MRI scanner and their MHD contaminated ECG was recorded and optically transmitted outside the faraday cage, then processed in real time using a previously developed software [START_REF] Abi-Abdallah | Cardiac and respiratory MRI gating using combined wavelet sub-band decomposition and adaptive filtering[END_REF].

B. Temporal and spectral analysis of the MHD contaminated ECG signals

For each subject, the one minute long signals recorded without any magnetic exposure and while exposed to a 1.5 Tesla static magnet were compared in both time and frequency domains. For temporal characterization we measured the amplitudes of the T waves and the R peaks, whereas in the frequency domain, the spectral density of each signal was estimated, and some characteristic values, were calculated to evaluate the spectral modifications induced by magnetic exposure.

C. Trigger production algorithm

Power spectral analysis of the ECG show that, P and T waves only present a significant spectral density up to 10 Hz, and that most of the QRS power lies in the 3-20 Hz band, beyond which the complex energy decreases gradually [START_REF] Thakor | Estimation of QRS Complex Power Spectra for Design of a QRS Filter[END_REF]. Thus, since gating relies on the R peaks, we can define this frequency band as being our region of interest in which we have all necessary information for trigger production. This concept was exploited in a previous work [START_REF] Abi Abdallah | Reference signal extraction from corrupted ECG using wavelet decomposition for MRI sequence triggering: application to small animals[END_REF] where we developed a peak detection algorithm based on wavelet sub-band decomposition. It consists in decomposing the ECG signal into multiple scales, then reconstructing the details and keeping only those that contain the maximum of the QRS energy in order to compose a reference signal. For example if the acquired signal had been sampled at 1Khz, we would decompose it into 8 scales and then reconstruct the reference signal by summing details d 6 and d 7 resulting in a [3.91-15.63] Hz sub-band. This reference signal is then subjected to a double threshold comparator to produce a cardiac trigger. The algorithm was tested on simulated as well as small rodents ECGs that had been contaminated by one of three imaging sequences (Gradient Echo, Fast Spin Echo and Inversion Recovery Spin Echo) and some analyzing wavelets revealed having outstanding performances according to the signal and the sequence type. The real time version of this algorithm uses adaptive filtering where a filter is computed to imitate the wavelet method performances with minimum delay. It has proven to be very efficient for insuring good cardiac image quality for small rodents [START_REF] Abi-Abdallah | Cardiac and respiratory MRI gating using combined wavelet sub-band decomposition and adaptive filtering[END_REF].

III. RESULTS AND DISCUSSION

A. Temporal and spectral characteristics

After analysis we found that the MHD artifact was most prominent on the lead I acquired signals; therefore in the following we only show numerical results of these signals. All recorded ECG signals showed a clear elevation of the T wave when the subject was placed inside the MRI scanner. Fig. 1 illustrates an example of two ECG signals of the same subject, acquired at 0 Tesla and 1.5 Tesla. For all subjects we calculated the T elevation as being the ratio of the T amplitude at 1.5 Tesla to the one measured at 0 Tesla, and evaluated the T/R amplitude ratios in each case (Fig. 2). The average T wave increase was found to be greater than 315%, where the T amplitude exceeds that of the R peak reaching over 1.2 times the R amplitude value. Note that, when the ECG is recorded without any magnetic exposure, the T wave amplitude is only around half that of the R peak.

For one subject the T wave amplitude even reached over twice that of the R peak, with approximately a 600% elevation. Fig. 3 shows an example of a spectral density estimate and Table 1 details the computed mean values and standard deviations of some characteristic spectral features. Globally we observed a spectral shift towards low frequencies, and an energy increase: the mean and median frequencies are reduced to almost half their values while the signal power doubles. On the other hand, the frequency where maximum PSD is attained remained relatively unchanged. It is noteworthy to mention that none of these magnetically induced changes were permanent. As stated by Gaffey [START_REF] Gaffey | Alterations in the rat electrocardiogram induced by stationary magnetic fields[END_REF] Tenforde [START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF] and Chakeres [START_REF] Chakeres | Effect of stati magnetic field exposure of up to 8 Tesla on sequential human vital sign measurements[END_REF] the ECGs recorded post exposure were similar to those recorded before exposure.

B. Evaluation of the trigger production algorithm

While the subjects were inside the MRI scanner tunnel their recorded ECG was processed using the real time trigger production algorithm described above in order to evaluate its performances in the presence of high MHD artifacts. When the algorithm was first ran while keeping the usual [3.91-15.63] Hz band to extract a reference signal used to produce the trigger, many false detections were made. Nevertheless if the lower band limit was pushed to 6Hz instead of 3Hz perfect detection would be attained. In fact, as can be seen in fig. 3 , beneath 5 or 6Hz the MHD induced spectral content would still be very dominant, thus in the extracted reference signal the T wave amplitude would remain high enough to be mistaken for an R peak, thus leading to false positive triggers. Fig. 5 illustrates an example of trigger production on a highly MHD contaminated ECG signal using the [6.64-26.56] Hz band to extract the reference signal. After processing with the sym7 wavelet all T waves were eliminated keeping only the R peaks for trigger production.

IV. CONCLUSION

The perturbation of the ECG signals collected during MR imaging greatly obstructs synchronization tasks. Moreover, the higher the requirements for good space-time or frequency resolution are, the more important the difficulties they induce. Some NMR generated artifacts may overlap the ECG spectrum, and are hard to filter out, such as the MHD induced voltages. These voltages alter the ECG waveform by elevating the amplitude of the T wave and might lead to R peak misdetections. In this work we examined the MHD induced artifacts on ECG signals recorded in a 1.5 Tesla magnet, and established a primary characterization of the induced changes using time and frequency domain analysis. T wave elevation rates were measured, and were found to agree with the values given by Tenforde. A great energy increase, concentrated in the low frequencies area, was noted. We also evaluated a previously developed real time algorithm for cardiac synchronization. The algorithm has been tested earlier on highly contaminated ECG signals for small rodents cardiac MRI gating, nevertheless the MHD effect was not very prominent in those signals given the relatively small size of rodents' arteries and low flow rate compared to those of humans. In this paper the method was tested on human ECGs with high MHD artifacts and has proven to be equally efficient by introducing a minor modification to the extraction process to narrow the extracted frequency band. This study will later be extended to include more subjects and greater magnetic field intensities to further investigate ECG alterations and test the efficiency of the reference extraction algorithm when even higher T wave amplitudes are involved. We will also conduct more experiments to test whether the modification of the reference signal band might affect the elimination of the gradient artifacts. Fig. 5. R peak detection on the contaminated ECG of subject n°5 performed while inside the MRI scanner using the sym7 wavelet. Despite the tremendous T wave elevation which amplitude reaches more than twice the value of the R peak (T/R=225%), the trigger production is 100% accurate when extracting a reference signal in the [6.64-26.56] Hz band.

Fig. 1 .

 1 Fig.1. Two ECG signals recorded on subject n°1, the first without any magnetic field influence and the second while exposed to a 1.5 Tesla magnet. We can clearly depict a prominent elevation (~280%) of the T wave which exceeds the R peak by approximately 20%

Fig. 4

 4 Fig. 4 plots the scalograms of the ECGs recorded on subject n°1 at 0 and 1.5 Tesla. It clearly shows a low frequency energy increase at the T waves locations.It is noteworthy to mention that none of these magnetically induced changes were permanent. As stated by Gaffey[START_REF] Gaffey | Alterations in the rat electrocardiogram induced by stationary magnetic fields[END_REF] Tenforde[START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF] and Chakeres[START_REF] Chakeres | Effect of stati magnetic field exposure of up to 8 Tesla on sequential human vital sign measurements[END_REF] the ECGs recorded post exposure were similar to those recorded before exposure.

Fig. 2 .

 2 Fig. 2. The 'T elevation' chart traces the ratios of the T wave amplitudes at 1.5 Tesla to those measured at 0 Tesla. The 'T/R 0T' and 'T/R 1.5T' charts show the ratios between the amplitudes of the T waves and those of the R peaks measured in the same signal when recorded at 0 or 1.5 Tesla respectively.

Fig. 3 .Fig. 4 .

 34 Fig.3. Power spectral density estimates of the ECG signals of subject n°1 recorded at 0 and 1.5 Tesla. One can observe a spectral shift towards low frequencies and a great increase of signal energy when the subject is exposed to a magnetic field.

TABLE I SPECTRAL VARIATIONS Mean Freq (Hz) Median Freq (Hz) Max Freq (Hz) Signal Power (V²)

 I

		0 T	1.5T	0 T	1.5T	0 T	1.5T	0 T	1.5T
	M	8.40	4.69	7.40	3.39	3.19	2.82	1.23	2.81
	SD	0.92	1.28	1.50	1.45	2.42	1.79	0.52	1.32

M : Mean, SD: Standard deviation