D Fraix-Burnet

E Davoust
email: davoust@ast.obs-mip.fr

C Charbonnel
email: corinne.charbonnel@unige.ch

The environment of formation as a second parameter for globular cluster classification

Keywords: Globular clusters, general -Galaxy, formation -Galaxy, evolution -Galaxy, halo -method, statistical

We perform an evolutionary multivariate analysis of a sample of 54 Galactic globular clusters with high-quality colour-magnitude diagrams and well-determined ages. The four parameters adopted for the analysis are: metallicity, age, maximum temperature on the horizontal branch, and absolute V magnitude. Our cladistic analysis breaks the sample into three novel groups. An a posteriori kinematical analysis puts groups 1 and 2 in the halo, and group 3 in the thick disc. The halo and disc clusters separately follow a luminosity-metallicity relation of much weaker slope than galaxies. This property is used to propose a new criterion for distinguishing halo and disc clusters. A comparison of the distinct properties of the two halo groups with those of Galactic halo field stars indicates that the clusters of group 1 originated in the inner halo, while those of group 2 formed in the outer halo of the Galaxy. The inner halo clusters were presumably initially the most massive one, which allowed the formation of more strongly helium-enriched second generation stars, thus explaining the presence of Cepheids and of very hot horizontal-branch stars exclusively in this group. We thus conclude that the "second parameter" is linked to the environment in which globular clusters form, the inner halo favouring the formation of the most massive clusters which subsequently become more strongly self-enriched than their counterparts of the galactic outer halo and disc.

INTRODUCTION

Globular clusters are touchstones of astrophysics.galaxies. The oldest of them witnessed the formation and early evolution of their host galaxies and of their substructures, and their study has historically coloured the different scenarios of galaxy formation. However, as a collective population in a galaxy, they present unsolved problems. In particular their origin is not firmly established despite the large amount of work devoted to the analysis of correlations among their observable properties. It has long been realised that part of the difficulty arises from sheer dynamical evolution undergone by these objects since the time of their formation. Indeed, any star cluster is the subject of a long list of erosive mechanisms that operate at different rates depending on the cluster's location and orbit within the Galaxy, and on its initial mass (Djorgovski & Meylan (1994), Gnedin & Ostriker (1997)).

Related to these difficulties, the search for a "second parameter", beyond metallicity, to explain the distribution of stars along the horizontal branch, has met with a limited success. [START_REF] Zinn | The Globular Cluster-Galaxy Connection[END_REF] suggested that the halo globular clusters break down into two groups according to their horizontal-branch properties: the two groups have different ages, kinematics and radial distributions. Rey et al. (2001) found that an age difference can explain different horizontal-branch morphologies at a given metallicity. But Dotter (2008), following the seminal work by Rood (1973) on the impact of stellar mass loss on the horizontal-branch morphology, has shown that considering α-element enhancement and metallicitydependent mass loss along the red-giant phase produces similar effects. On the other hand Lee, Gim & Casetti-Dinescu (2007) have found that globular clusters with extended horizontal branch are more massive than normal clusters and are dominated by random motions with no correlation between kinematics and metallicity. Multivariate analyses have been performed by Fusi Pecci et al. (1993) who find that more concentrated clusters have bluer and longer horizontal-branch, and by Recio-Blanco et al. (2006) who find that more massive clusters have a horizontal branch that extends to higher temperatures. ronment at birth must influence other cluster parameters. This was already recognized by Zinn (1985), who combined metallicity and kinematics to establish two major groups (halo versus bulge/disk). The properties of the horizontal branch have been used, but its usual characterisation by the parameter HBR = (B-R)/(B+V+R) is also unsatisfactory. In this respect, the parameter T e (maximum effective temperature on the horizontal branch) introduced by Recio-Blanco et al. ( 2006) is a welcome innnovation as shown in the present paper.

The classical paradigm describing globular clusters as fairly simple systems of coeval stars of homogeneous chemical composition has been seriously challenged recently, and this may bring a crucial piece to the puzzle. One fundamental characteristic of these systems is their metallicity [Fe/H] (identified as the so-called "first-parameter" (van den Bergh (1967))) that is generally inferred from their integrated photometric colours, and that varies strongly from cluster to cluster; in our Galaxy, globular clusters have [Fe/H] ranging between ∼ -2.2 and 0 (Harris (1996)). Spectroscopy reveals that, within individual clusters, stars present very homogeneous contents in Fe, but also in α-and s-elements, indicating that protoclusters formed from gas pre-enriched in heavy metals (James et al. (2004); [START_REF] Prantzos | [END_REF]). This is in agreement with the predictions of quantitative models that rule out stochastic self-enrichment in most globular clusters as a significant contributor to their heavy metals, leaving pre-enrichment as the dominant contributor to metallicity [Fe/H] ( Bailin and Harris (2009)).

However globular cluster stars exhibit extremely scattered light-element (Li, F, C, N, O, Na, Mg, and Al) contents that are not seen among their field counterparts (see e.g. Carretta (2006)). This points to early internal chemical evolution (i.e., self-enrichment) in the globular cluster driven by first-generation massive and fast evolving stars which polluted with their hydrogen-burning products (among which helium in very important quantities) the intracluster gas out of which second generation stars formed (Decressin, Charbonnel & Meynet (2007), [START_REF] Prantzos | [END_REF]). Recent findings of double or even multiple stellar populations in the colour-magnitude diagrams of several globular clusters, as well as the complexity of the horizontal branch morphology (namely the wide colour distribution, i.e., effective temperature, of the stars presently burning helium in their core) constitute further evidence for internal evolution [START_REF] Piotto | IAU Symp[END_REF]). All these features can indeed be related to the presence of a second generation of He-enriched stars. Importantly, the star formation history depicted by these features seems to vary from cluster to cluster (Milone et al. (2008)) in a way which is still far from being understood. However we have now firm evidence that Galactic globular clusters have undergone internal chemical evolution and complex star formation histories during their infancy that shaped their properties and in particular their present total mass (Decressin, Charbonnel & Meynet (2007)). This new paradigm has opened a novel route for a better understanding of the origin and history of globular clusters.

It thus appeared to us that a multivariate analysis which simultaneously takes into account any cosmic variance due to evolving physical conditions and groups objects according to environment of formation would be very valuable. Cladistics provides such a methodology. It differs from other clustering analyses in that it focuses on evolution within and between groups rather than on similarities between objects (Wiley et al. (1997)). Cladistics is very commonly used in evolutionary biology and has been pioneered in astrophysics by Fraix-Burnet et al. (2006a,b); [START_REF] Thuillard | [END_REF] and successfully applied to the dwarf galaxies of the local group (Fraix-Burnet et al. (2006c)). This paper presents a multivariate analysis based on the method of cladistics of a large sample of Galactic globular clusters. After presenting the data and the method of analysis (Sect. 2), we describe the three groups found by the cladistic analysis (Sect. 3.1), and discuss two important results, evidence for self-enrichment (Sect. 4) and a possible luminosity-metallicity relation (Sect. 5). We then compare the properties of the three groups with those of Galactic halo field stars (Sect. 6), before proposing a scenario for the formation of the three groups (Sect. 7).

DATA AND METHOD OF ANALYSIS

The choice of parameters is a crucial step in any multivariate analysis. Djorgovski & Meylan (1994) have shown that the manifold of Galactic globular cluster properties has a dimension larger than 4, but that a subset of parameters linked to morphology and dynamics forms a three-dimensional family. Including properties of the stellar populations (e.g. a horizontal-branch parameter or metallicity) will increase the manifold by 1 or 2 dimensions (Fusi Pecci et al. (1993)). Finally, using a large number of photometric and structural parameters, Recio-Blanco et al. (2006) found that 4 eigenvectors account for 79% of the total sample variance.

Taking advantage of this indication, we selected the following four parameters for analysis: relative ages, metallicity ([Fe/H]), absolute V magnitude (M v ), and maximum effective temperature (T e ) on the horizontal branch. The age parameter is related to the secular evolution of the stellar populations. [Fe/H] reflects the chemical composition of the environment when and where globular clusters formed and is the "first parameter" for the horizontal branch morphology M v is a structural parameter that measures the present total baryonic mass of the globular clusters1 . Finally T e is a measure of both the pristine chemical composition of the protocluster ([Fe/H] being the first parameter) and of the helium enrichment during early internal chemical evolution, since stars with higher helium content are expected to reach higher effective temperatures on the horizontal branch (D'Antona et al. (2002); Recio-Blanco et al. (2006)). The last three quantities describe truly intrinsic properties of globular clusters. As the age parameter evolves in all clusters, it cannot be used to classify them at the same level as the other three parameters. We thus gave it a lower weight in the cladistic analysis (see Appendix A).

We performed our analysis using the large sample (54 objects) of Recio-Blanco et al. (2006) based on homogeneous Hubble Space Telescope photometry. We used the T e values obtained uniformly from this database by Recio-Blanco et al. (2006), as well as the relative ages and M v values they adopted (i.e., taken respectively from De Angeli et al. (2005) and the 2003 on-line revision of Harris (1996)). We did not include more parameters in the cladistic analysis as we preferred to avoid the unwanted effect of redundancies, which give more weight to correlated parameters. Note also that we did not use any kinematical information in the cladistic analysis. However we used other parameters a posteriori to characterize the different groups found by the cladistic analysis: the radial velocities and structural parameters were taken from the 2003 on-line revision of Harris (1996), the orbital parameters from Dinescu, Girard & van Altena (1999); Dinescu et al. (2003); Casetti-Dinescu et

The environment of formation as a second parameter 3 al. (2007). In this respect, the distinct orbital properties of group 3 found a posteriori (see §3) are independent of the methodology.

More recent age estimation was published by Marin-Franch et al. (2009) after most of this project was completed, but for only 35 out of 54 globular clusters of our sample. Using these produces an inhomogenous data set, from different sources, relying on different values of Fe/H. The ages of Marin-Franch et al. (2009) have been determined from colour-magnitude diagrams and values of Fe/H which are different from those of Recio-Blanco et al. (2006). Marin-Franch et al. (2009) themselves point out the importance of using a homogeneous set of Fe/H to derive the ages. Nevertheless, it is instructive to perform analyses using both sets and compare the results, so as to determine how sensitive they are to the specific choice of parameter values. One then has the problem of combining two sets of ages, and it is not obvious how this should be done. It turns out that 9 out of 11 globular clusters that calibrate relative ages of Marin-Franch et al. (2009) are in common with De Angeli et al. (2005), and using the relative ages of De Angeli et al. (2005) for these 9 globular clusters gives a mean value of 1.0055 instead of 1.00, which is fine. However, comparing the ages of all the globular clusters in common suggests a non-linear systematic effect, which should perhaps be taken into account. Furthermore, we cannot simply convert each set of relative ages to absolute ages with the zero-point of each set, because the two zero-points are rather different (11.2Gyr for De Angeli et al. (2005) and 12.8Gyr for Marin-Franch et al. (2009)). We thus simply used the relative ages without any attempt at homogenizing them, and the zero-point of one author for all ages. This additional analysis is compared to the main one in Sect. 3.2.

The multivariate analysis was performed using the method of cladistics. In short, the method works as follows. One first builds a matrix with values of the four parameters for each clusters. The values must be discretized, and the number of bins (here 10) depends on the resolution one wants for the analysis. One then chooses a cluster which represents the most unevolved state, in the present case the metal-poorest cluster (NGC 6934), and the software classifies all the other clusters in order of increasing diversification of properties (in other words, by increasing distance in the manifold of parameters). Clusters which are diverging from the original cluster in the same direction are put on the same branch. We refer to Appendix A and Fraix-Burnet et al. (2006a,b,c);[START_REF] Thuillard | [END_REF] for more details on the principles of the method.

THREE GROUPS OF CLUSTERS

The main tree

The main result of our analysis is presented in the form of a tree structure, a usual form of representation in graph theory. The properties of the sample can be read from the structure of the tree shown on Fig. 1.

The tree has been rooted with group 2 which has the lowest metallicity on average, and as such is supposedly made of more primitive (or "ancestral") material. The tree divides into three main branches, which define three groups with quite distinct properties. Age increases roughly monotonically along each branch, as expected. There are some subbranches sharing similar values of the adopted parameters, such as NGC 5904, 2808, 6388 and 6441, and which set them apart in the four-dimensional space of parameters, implying that they might in fact belong to a fourth small group with properties similar to those of group 1. Specialists will immediately Figure 1. The tree resulting from the cladistic analysis. The sample break into three distinct branches. Group 1 is composed of inner halo clusters, group 2 of outer halo clusters, and group 3 of disc clusters remark that the latter three clusters indeed share peculiar properties, in particular a very Helium-rich stellar population (Pumo, D'Antona & Ventura (2008)). We emphasize that the Helium abundance is not one of the parameters included in the cladistic analysis: it thus must influence one way or another the four parameters used in the analysis.

The properties of the three groups are presented in Table 1. The first 7 rows give the characteristics of the 3 groups: number of clusters in each group, mean distance from the Galactic center R gc , height above the Galactic plane Z, metallicity [Fe/H], absolute magnitude (M v ), mean velocity of rotation V rot in the Galactic plane (computed with the equations given by Frenk & White (1980) using a velocity of the Sun of 220 km/s), radial velocity dispersion σ. The next 6 rows describe the correlations between the four parameters. The bottom part of the Table lists orbital parameters taken from Dinescu, Girard & van Altena (1999); Dinescu et al. (2003); Casetti- Dinescu et al. (2007), which were not available for all clusters of each group: successively number of clusters, period of rotation P, total energy E, eccentricity of the orbit e, apocentric distance R a , maximum distance reached above the Galactic plane Z max , inclination angle with respect to the Galactic plane Ψ, angular momentum L, and finally the velocity components in cylindrical coordinates: vertical velocity |W|, radial velocity Π and tangential velocity Θ.

Paired t-tests showed that the differences of the means of the groups taken two by two is not equal to 0 (p < 0.05) for the parameters log(T e ), R gc , [Fe/H], e. This is also the case between group 1 and group 2 and between group 1 and group 3 for M v and Θ, and between group 2 and group 3 for Z and age. There is also evidence (0.05 < p < 0.1) for different means of the latter two groups for P, E, R a , Z max and |W|, as well as between group 1 and group 2 for age and between group 1 and group 3 for |W|. We emphasize that the rotational properties of Galactic globular clusters are very uncertain, since they are derived from projected radial velocities or numerical simulations of orbits in a model Galaxy. These properties show that the first two groups belong to the halo population of clusters, while the third group belongs to the thick disc population. Hereafter, the thick-disc clusters will simply be called disc clusters. The average velocity of rotation of group 1 and 2 together is V rot = 9 km/s. Group 3 is confined to the Galactic plane, and has a high V rot and low σ. If we separate group 3 into two subgroups of equal size according to their distance from the Galactic center, we find that V rot is 88 km/s for the inner subgroup (R gc < 3 kpc) and 187 km/s for the outer subgroup. There is also evidence that group 3 has a shorter P, lower e, Ψ and Z max , and no radial motion, as expected from clusters that partake in the overall rotation of the disc.

One cluster has certainly been misclassified. NGC 6981 is in group 3 although it is at Z = 9.1 kpc and has a low velocity of rotation. Since NGC 6981 is a borderline cluster in all the Figures, the value of one of the four parameters may be erroneous. Indeed, raising M v from -7.04 to -7.27 brings it into the next bin in our cladistic analysis (which requires that the data be discretized into a limited number of bins), and running the cladistic analysis again moves the cluster to group 2. Moving any of the other three parameters by one bin and redoing the cladistic analysis does not change the status of the cluster. Since we found no reason for an erroneous M v , we left it in group 3.

Another possible discrepant cluster is NGC 6266, which is in group 1, but which, according to Dinescu et al. (2003), belongs to a rotationally supported system, on the basis of its kinematics (but without precise orbital determination), despite its low metallicity. However, it is on the same subbranch as (and undistinguishable from) NGC 7089, which definitely belongs to the halo, according to our analysis and that of Dinescu, Girard & van Altena (1999). In addition, as pointed out by the referee, an isotropic distribution of orbits will statistically produce one or several ones in or near the Galactic plane. We are thus confident that NGC 6266 belongs to the halo.

We now compare the statistical properties of the two halo groups with those of the Galactic halo field stars. The dichotomy of the Galactic halo stellar population has been suspected for some time. The most quantitative study in that respect, that of Carollo et al. (2007), clearly identifies two broadly overlapping structural components corresponding to an inner and an outer halo. Stars of the inner halo are in highly eccentric orbits, in slightly prograde rotation, and have an average metallicity of [Fe/H] = -1.6. The outer halo stars have a uniform distribution of eccentricities, are in highly retrograde orbits, and of lower metallicity [Fe/H] = -2.2. These properties are also among those that distinguish the two groups of halo clusters, and indicate that the clusters of group 1 may have originated in the inner halo: they have higher eccentricities and metallicities than group 2, while group 2 formed in the outer halo: they have the largest R gc , R a , Z and Z max , lowest metallicities; Θ is positive, but not significantly so.

Additional analysis using an inhomogeneous set of ages

We also applied our analysis to the sample of 54 globular clusters with composite ages as explained in Sect. 2. It leads to roughly the same three groups as before, with several obvious halo clusters moving into G3 (which is in principle composed of thick-disk globular clusters), and two globular clusters of G1 moving into G2, among them NGC 6218 which has a Cepheid (and should thus be in G1). In other words, the number of misclassified globular clusters rises from one (NGC 6981) to only about 6. This indicates that our groups are fairly robust. Of course, their detailed contours depend on the choice of data and we believe that the slight discrepancy is due to the inhomogeneous data of this additional analysis: the ages of Marin-Franch et al. (2009) have been determined from colour-magnitude diagrams and values of Fe/H which are different from those of Recio-Blanco et al. (2006), while logTe comes from Recio-Blanco et al. (2006) and is derived from diagrams using the corresponding set of Fe/H values.

EVIDENCE IN FAVOUR OF SELF-ENRICHMENT

Self-enrichment by a first generation of stars is frequently advocated to explain the chemical anomalies in globular clusters (e.g. [START_REF] Prantzos | [END_REF]) in relation with horizontal-branch morphology (D'Antona & Caloi (2008)). We list below distinctive properties of the three groups of globular clusters pointing to such a process.

The well-known correlation between metallicity and extent of the horizontal branch is clearly present in Fig. 2: the more metallic the cluster, the less extended its horizontal branch at a given metallicity. Also, the more luminous (that is the more massive) clusters tend to have more extended horizontal branch (Recio-Blanco et al. (2006); Lee, Gim & Casetti-Dinescu (2007)). This latter point is easily understood in the self-enrichment framework, in the sense that more massive globular clusters retain the helium-rich ejecta of massive polluter stars in their deeper potential well more efficiently than less massive globular clusters .

The crucial new information brought by our analysis in this context is that the present mass of the clusters seems to depend on their origin, the inner halo clusters being presently more massive than their outer halo counterparts (Fig. 2). In fact, in order to have extreme light-element abundance patterns, and in particular extreme O-Na anticorrelation (Carretta (2006)), which is linked to the extent of the horizontal branch (Carretta et al. (2007)), the inner halo clusters must have been even more massive in the past, before they lost a huge number (96%) of first-generation low-mass stars in their early dynamical evolution (Decressin, Charbonnel & Meynet (2007)). We thus expect the mass difference between the inner and outer halo globular clusters to have been even larger in the past.

The disc and halo clusters are well separated in the agemetallicity diagram, shown in Fig. 3. For the halo clusters, metallicity decreases with age. For the disc clusters, on the contrary, it marginally increases with age, if at all. The spread in age is 4 Gyr for the halo component and only 1.5 Gyr for the disc one. The figure confirms that the metallicity of NGC 2808 should indeed be about 0.5dex lower. The two other He-rich clusters of group 1 do not appear on this plot because no age estimate is available for them, but their metallicity is indeed about 0.5dex higher than the highest metallicity of the rest of the halo population. This supports the claim by Caloi & D'Antona (2008) and [START_REF] Prantzos | [END_REF] (see also Decressin, Charbonnel & Meynet (2007)) that He-enrichment must be associated with the build-up of abundance anomalies of light elements during the phase of self-enrichment.

The presence of multiple stellar populations in the colourmagnitude diagram is another evidence for self-enrichment. Unfortunately, only three of the clusters of the sample, NGC 1851, 2808 and 6388, all belonging to group 1, are known to have such a feature [START_REF] Piotto | IAU Symp[END_REF]). We predict that two other known such clusters, NGC 5139, 6656 which are both metal-poor and massive, are also inner halo globular clusters.

Additional clues to the self-enrichment scenario can be gathered from the RR Lyrae and Cepheid contents of the three groups. Globular clusters have historically been divided into two groups (Oosterhof I and II) according to the properties of their RR Lyrae stars. Using the compilation of Clement et al. (2001), we find that the two Oosterhof types are equally present in group 1 and 2. More interestingly, we find that the distribution of periods of RR Lyrae stars in group 2 is more sharply peaked than the corresponding distribution for group 1, which presents a minor secondary peak at a higher period. The narrower period distribution in group 2 implies a small dispersion in mass loss along the red-giant branch (Caloi & D'Antona (2008)), while the wider distribution in group 1 implies several generations of stars, each with a narrow distribution of mass loss along the red-giant branch, and with increasing helium abundance (D'Antona & Caloi (2008)), reinforcing the necessity of self-enrichment.

Furthermore all the population II Cepheids are found in clusters of group 1, and none in clusters of group 2. Population II Cepheids result from the evolution of post-horizontal branch stars which start from the higher temperatures of the zero-age horizontal branch and move toward the asymptotic giant branch or leave that branch on rapid blueward loops (Wallerstein (2002)); this explains their absence in halo clusters with low T e.

We thus reach the conclusion that the inner halo favours the formation of very massive clusters, which retain more easily the products of first-generation stars and thus become more strongly self-enriched, giving rise to more extended horizontal branches.

A LUMINOSITY-METALLICITY RELATION FOR GLOBULAR CLUSTERS?

A mass-metallicity or luminosity-metallicity relation is found among galaxies, but is not expected and has not been found in the Galactic globular cluster system (Djorgovski & Meylan (1994)), nor in numerical simulations of globular cluster formation (Kravtsov & Gnedin (2005)). However, it has been found, in the form of a "blue tilt", in colour-magnitude diagrams of globular cluster systems in bright ellipticals (Brodie & Strader (2006); Harris et al. (2006); Harris (2009); Mieske et al. (2006)). Such a relation, where metallicity increases with Mv, is present in our sample (Fig. 4), if we consider the disc and halo clusters separately.

We have plotted a line separating the disc and halo clusters in Fig. 4, which can be used as a criterion for distinguishing the two types of clusters, in conjunction with other criteria, since the separation is not perfect (NGC 104 is a notable exception). The disc clusters are on average fainter than the halo ones by about 3 magnitudes at a given metallicity. Correlation lines were adjusted to the two subsamples: they have comparable slopes of -2.0 and -2.8 for the halo and disc clusters respectively, with admittedly low correlation coefficients of 0.53 and 0.32. These slopes are much lower than that of -6.75 predicted and found for the dwarf galaxies of the Local Group (Dekel & Silk (1986)). They were however probably much steeper when the globular clusters formed, since we expect some of these objects to have lost a large fraction of their initial mass in the self-enrichment framework.

A remarkable property of our disc clusters is that they extend to metallicities lower than the conventional limit of [Fe/H] = -0.8 (Zinn (1985)). This is not very surprising per se, since Dinescu, Girard & van Altena (1999) have found three metal-poor Galactic globular clusters with thick-disc kinematics, and in M31 there is also evidence for metal-poor globular clusters with disc kinematics (Morrison et al. (2004), although see Fusi Pecci et al. (2005)).

If we divide group 3 along the conventional limit, we find that the metal-rich clusters have V rot = 184 km/s whereas the other ones have V rot = 71 km/s, rather low, but still significantly larger than that of the halo clusters. We have checked that this low mean velocity is not due to one cluster in particular. The two subgroups do not distinguish themselves otherwise; in particular they have the same spatial distribution (same mean R gc and Z). Two of the low-metallicity disc clusters, NGC 6171 and 6362, have orbits determined by Dinescu, Girard & van Altena (1999), which confirm that they do belong to the disc population. In fact Dinescu, Girard & van Altena (1999) state that their most significant result is to have shown the existence of metal-poor clusters with orbits consistent with the thick-disc motion. Our analysis confirms this finding.

The second important result of the present paper is that the disc and halo globular clusters should not be separated on the basis of metallicity, but rather of a multivariate analysis, taking into account other parameters. We propose to use the magnitude at a given metallicity as a rough criterion, with a limit such that M v = -2.4 × [Fe/H] -10.45, together with other criteria, such as location in the Galaxy and velocity of rotation.

COMPARISON WITH OTHER STUDIES AND WITH HALO FIELD STARS

Before interpreting the differences between the three groups in terms of formation history, we compare them to the traditional disc, young halo and old halo groups of [START_REF] Zinn | The Globular Cluster-Galaxy Connection[END_REF], to emphasize that they are rather different. There are 4, 4 and 1 young halo clusters in group 1, 2 and 3 respectively, and the numbers are 17, 10 and 8 for the old halo clusters. Two of Zinn's disc clusters (NGC 6388 and 6441) are in our group 1. Since they are located near the Galactic centre and have no determined orbits, it is not possible to decide if they kinematically belong to the disc or the halo. Furthermore, as mentioned in Sect. 3.1, these two clusters might in fact belong to a fourth small group, maybe of Galactic bulge clusters. Turning to the metallicity vs HBR diagram, we find that group 1 extends to lower metallicities at a given HBR than the old halo clusters, and that there are clusters of group 1 and 2 among the old halo clusters with the reddest HBR. Comparing our groups to those of Lee, Gim & Casetti-Dinescu (2007), we find that all the clusters of their group with extended horizontal branch are in our group 1, except NGC 4833 (which is in our group 2).

If we now compare (Table 2) our grouping with that of Harris (2001) (see his Table 1.6), we find that G3 dominates in metal-rich clusters class (MRC) and G1 and G2 dominate in metal-poor cluster class (MPC). In MPC alone there is no clear separation in R gc between G1 and G2, while G3 tends to be in the inner regions. G2 tends to be among the more metal-poor globular clusters and tends to be at larger galactocentric distances. In summary, the MPC/MRC dichotomy corresponds roughly to our halo/disk separation, and our G2 populates the very metal-poor and distant MPC. Table 2 also shows that a classification based on arbitrary criteria does not quite retrieve the groups obtained with a multivariate analysis.

We now compare the content in α-elements of the two groups of halo globular clusters with those of halo field stars. Several studies have shown that the field stars in the outer halo, identified through their kinematical or orbital parameters, tend to have lower and/or more dispersed relative abundances in α-elements (Stephens & Boesgaard (2002); Gratton et al. (2003)), which points to a difference in star formation rates of their birth environment: the αelements are indeed almost exclusively provided by core-collapse supernovae, which arise from ephemeral massive stars, while Fe is essentially produced in type Ia supernovae, which arise from stars with longer lifetimes. We find marginal evidence for a similar difference between groups 1 and 2 (see Fig. 5). [α/Fe] decreases with increasing [Fe/H] in group 1, while it is more dispersed and shows no clear trend with [Fe/H] in group 2. This confirms a similar origin for the α-elements in globular clusters and field stars of the halo, which thus has to occur prior to protocluster formation.

ORIGIN OF THE THREE GROUPS OF GLOBULAR CLUSTERS

It is generally assumed that the Galaxy assembled in a hierarchical fashion from collapsing haloes of dark matter (Bertschinger (1998).

Small proto-galactic clumps formed first, from initial small-scale density fluctuations, and collapsed in a dissipationless way (because their gas was quickly consumed or blown away), or else they merged and grew in size to form larger clumps which spiralled toward the inner regions of the Galaxy by dynamical friction and experienced a dissipational collapse. This scenario, or variants of it, first proposed by Searle & Zinn (1978) invoked to explain the properties of the Galactic globular clusters and of the stellar halo. Furthermore, cosmological simulations have shown that globular clusters can form in the densest regions of collapsing subhaloes of dark matter and giant molecular clouds, when the clouds reach a critical density and are under a high external pressure. The mass distribution function of the clusters is similar to that of the clouds (Kravtsov & Gnedin (2005)).

It has often been proposed that some clusters, in particular those identified as young halo clusters by [START_REF] Zinn | The Globular Cluster-Galaxy Connection[END_REF], which tend to be in the outer halo and counter-rotating, were formed by accretion and disruption of satellite galaxies. But the chemical homogeneity of the halo, as well as substantial differences in chemical composition between field stars in the halo and dwarf spheroidal galaxies, argue against the accretion scenario (Stephens & Boesgaard (2002); Pritzl, Venn & Irwin (2005);Geisler et al. (2007)).

The properties of our three groups of clusters can be interpreted in the following way, without resorting to an external origin for any of the groups:

-The clusters of the outer halo (group 2) formed during the initial dissipationless collapse of the proto-galaxy, from material already polluted by earlier generations of stars, but not well homogenized and thus inhomogeneous in α-elements. Contrary to the outer halo stars, they lost their initial average retrograde rotation by dynamical friction and gravitational encounters. As suggested by the referee, this group could also have originated in some "pregalactic dwarfs" (i.e. metal-poor, gas-rich satellites that soon afterward began hierarchical merging).

-The clusters of the inner halo (group 1) formed later, during the dissipational phase of Galactic collapse, which continued in the halo after the formation of the thick disc and its globular clusters. Since the formation of group 1 occurred later, the molecular clouds from which they formed had time to grow by accretion of smaller clumps. These clouds were already enriched at the same level in α-elements. Thanks to the strong potential well in the clusters (as evidenced by their high central velocity dispersions), the He-rich ejecta of first generation massive stars were not blown away and found their way into a more strongly helium-enriched second generation of stars, favouring the production of hot horizontal-branch and Cepheid stars.

-As indicated by their short range in age (1.5 Gyr, see Fig. 3), the disc clusters (group 3) formed in a more rapid fashion than the two other groups, before many clusters of group 1. This could presumably be due to the higher densities and external pressure in the thick disc. This group shows significant average prograde rotation, because the dissipational collapse of the disc conserved angular momentum. The metal-poor disc clusters seem to rotate more slowly and have larger eccentricities and inclinations than the metal-rich disc clusters. Since there is no significant age difference between the two subgroups, we assume that the metal-poorer clusters formed further away from the Galactic plane, and thus retained a larger vertical velocity component. use of cladistics is justified a priori by the evolutionary nature of globular clusters, and a posteriori by the strong astrophysical significance of the grouping found.

In astrophysics, cladistics has already been applied to galaxies because they can be shown to follow a transmission with modification process when there are transformed through assembling, internal evolution, interaction, merger or stripping (Fraix-Burnet et al. (2006a,b,c);[START_REF] Thuillard | [END_REF]). For each transformation event, stars, gas and dust are transmitted to the new object with some modification of their properties. For globular clusters, interactions and mergers are probably rare. It was previously thought that once they assembled, only the stellar ageing would affect their properties. Nowadays, we have firm evidence that internal evolution can create another generation of stars, and clusters can lose mass. Basically, the properties of a globular cluster strongly depend on the environment in which it formed (chemical composition and dynamics), and also on the internal evolution which includes at least the ageing of its stellar populations. To compare globular clusters, it is thus necessary to take into account the different stages of evolution of both the objects and their environments of formation. Since the clusters form in a very evolving environment (evolution of the Universe and the dynamical environment of the parent galaxy), the basic properties of different clusters are related to each other by some evolutionary pattern. In particular, the dust and gas from which the stars of the globular clusters form have been "polluted" (enriched in heavy elements) by more ancient stars, being field stars or belonging to other globular clusters. This results in a kind of transmission with modification process, which justifies a priori the use of cladistics. It must be clear that this is not a "descent with modification" in the sense that there is no replication. But evolution does nevertheless create diversity. We are dealing with phylogeny (relationships between species), not with genealogy (relationships between individuals). Since a multivariate classification of globular clusters is not yet available, we assume in the present work that each cluster represents a species that will have to be defined later on.

As our work on galaxies has shown us, it is important to remove parameters that are redundant. Since previous studies of the manifold of Galactic globular clusters have shown that 4 parameters are sufficient to describe their diversity, we selected 4 parameters, 3 of which are intrinsic characteristics of the environment of formation. The fourth one, age, is particular in the sense that it does not inform on the conditions when the clusters formed, and is not discriminant for clustering because it evolves similarly for any cluster (parallel evolution). However, age is useful to rank the clusters within each group. Consequently, we applied to age a weight half that of the other three parameters. In addition, a stepmatrix was employed to impose the irreversibility constraint on the age parameter (age can only increase). In contrast to multivariate distance methods, undocumented values are not a problem in cladistics analyses. This is why the seven galaxies that have no age determination (see Fig. 3 and Table 1) have not been excluded in our work.

In this paper, we use parsimony as the optimisation criterion. This works as follows. One first builds a matrix with values of the four parameters for all clusters. The values for each parameter are discretized into 10 bins representing supposedly evolutionary states. Discretization of continuous variables is quite a complex problem, especially in the evolutionary context (see e.g. Goloboff et al. (2006); [START_REF] Thuillard | [END_REF]). The choice of the number of bins cannot be made in a simple objective way. Here, we took equal-width bins, and considered a compromise between an adequate sampling of continuous variables and the uncertainties on the measurements. The first constraint is given by the software (32 in this case). The second one would a priori give a lower limit of something like total range / uncertainty, but Shannon's theorem would multiply this by 2. Hence, 10 bins would account for about 20% measurement errors, which is rather large, but Recio-Blanco et al. (2006) do not provide precise error estimates, especially for logTe. Even so, border effects always imply that some objects could belong to a bin or its neighbour, a process that add some more artificial noise. The best way to avoid this effect is to make several analyses with different number of bins and check that the result does not depend on this number. We have done this for 3, 5, 8, 10, 12, 15 and 20, and the result is identical, to at most one misplaced cluster, for numbers higher than 8. For 3, no structure is found, and for 5 bins the groups are not well defined.

Then, all possible arrangements of clusters on a tree are constructed, and using the discretized matrix, the total number of state changes is computed for each tree. The most parsimonious tree is finally selected. If several such trees are found, then a consensus (strict or majority rule) tree is built. The whole procedure is computerised since the number of arrangements is here very large. The result is a diversification scenario that should be confronted to other knowledge and parameters. Maximum parsimony heuristic searches were performed using the PAUP*4.0b10 (Swofforf ( 2003)) package. The results were interpreted with the help of the Mesquite software [START_REF] Maddison | Mesquite: a modular system for evolutionary analysis[END_REF]).

The tree presented in Fig. 1 is a majority rule consensus tree of 20000 trees, the strict consensus tree showing exactly the same three groups but with group 1 and 2 slightly less resolved. To further assess the robustness of the tree, it was not possible to make bootstrapping due to the irreversibility constraint on the age parameter, and it would not have been very significant with only 4 parameters. We performed other analyses using only 3 parameters, excluding the age. They all gave essentially the same three groups, but they were individually slightly less resolved, as expected. All these convergences yield strong confidence on the tree shown in Fig. 1. In the end, the most important point is the astrophysical interpretation we are able to give of the results.

On Fig. 1, the tree is rooted with group 2 as an outgroup. This is not strictly necessary in cladistics, and we find here the same three groups whatever the root chosen or even on the unrooted tree. But we know that a low metallicity is an ancestral state for stars in general, this is why we have chosen group 2 a posteriori because it has a homogeneously low metallicity.

Figure 2 .Figure 3 .

 23 Figure 2. Metallicity-logTe diagram with symbol size indicating visual magnitude. The halo clusters are represented by open symbols: squares for group 1 and triangles for group 2. The disc clusters (group 3) are represented by full squares

Figure 4 .

 4 Figure 4. Metallicity-M v diagram. Same symbols as Fig. 2. The solid line is the limit between disc and halo clusters according to our criterion (M v = -2.4×Fe/H -10.45). The disc cluster above this line is NGC 104

Table 1 .

 1 Properties of the three groups of globular clusters. No age estimate is available for NGC 6139, 6229, 6304, 6388, 6441, 6569 and 6642. The middle part of the Table gives the correlations: + means a correlation,means an anticorrelation, x no correlation. The orbital properties presented in the lower part of the Table are only available for a subset of each group.

	Numbers in brackets are rms dispersions
		Group 1	Group 2	Group 3
		inner halo	outer halo	thick disc
	Number of clusters	25	11	18
	R gc (kpc)	9.4 (7.4)	12.9 (8.0)	4.2 (2.9)
	Z (kpc)	4.8 (4.5)	8.6 (7.6)	1.9 (2.0)
	Fe/H	-1.40 (0.35)	-1.92 (0.16)	-0.92 (0.35)
	M v (mag)	-8.5 (0.7)	-7.6 (0.6)	-7.1 (0.9)
	V rot (km/s)	-7.	+46.	+119.
	σ (km/s)	120 (107)	151 (107)	69 (74)
	Age	9.98 (0.96) a	11.17 (0.70)	10.18 (0.48) b
	Age-log(T e )	+	x	-
	Age-[Fe/H]	-	-	+
	Age-M v	x	x	x
	log(T e ) -[Fe/H]	+	x	+
	M v -[Fe/H]	+	+	+
	M v -log(T e )	+	x	x
	Number of clusters	12	8	5
	P (Myr)	353 (212)	391 (254)	142 (35)
	E (10 2 km 2 /s 2 )	-691 (341)	-649 (339)	-1027 (111)
	e	0.63 (0.18)	0.54 (0.18)	0.21 (0.10)
	R a (kpc)	17 (12)	19 (13)	6 (1.3)
	Z max (kpc)	7.2 (6.2)	9.6 (8.5)	1.5 (1.0)
	Ψ (deg)	32 (12)	38 (17)	21 (15)
	L	886 (802)	941 (575)	866 (314)
	|W| (km/s)	100 (78) c	77 (40)	34 (16) d
	Π (km/s)	+19 (141)	-33 (124)	-3 (22)
	Θ (km/s)	+15 (121)	+89 (129)	+170 (33)

a average of 21 values; b average of 15 values; c average of 13 values; d average of 6 values.

  Same symbols as Fig. 2. No age estimate is available for NGC 6139, 6229, 6304, 6388, 6441, 6569 and 6642. Metallicity decreases with age in the halo clusters while it very marginally increases with age in the disc clusters

Table 2 .

 2 Metallicity versus light-element abundance. Same symbols as previous figures. The line is a least-squares fit to the globular clusters of group 1. Small triangles are disc (green) and halo (blue) field stars. Small red open squares are dwarf galaxies of the Local Group (fromVenn et al. 2004) Comparison between Table1.6 of[START_REF] Harris | Saas-Fee Advanced Course[END_REF] and our grouping by number of clusters for each class.

	and verified by cos-

The absolute magnitude M v of globular clusters in the Milky Way spans a vast range (-1.7 M v -10.2, Harris (1996)), and reflects a large mass range (10 3 -10 6 M ⊙ , McLaughlin & van der Marel (2005)). c

RAS, MNRAS000, 1-9

c 2009 RAS, MNRAS 000, 1-9

c 2009 RAS, MNRAS 000, 1-9 6 D. Fraix-Burnet, E. Davoust & C. Charbonnel

ACKNOWLEDGEMENTS

Paired t-tests were performed with the VOStat software (http://vo.iucaa.ernet.in/ voi/VOStat.html). We thank Asis Chattopadhyay for advice. We would like to thank the referee, W. Harris, for his constructive comments and suggestions.

APPENDIX A: CLADISTICS APPLIED TO GLOBULAR CLUSTERS

Multivariate clustering methods compare objects with a given measure and then gather them according to a proximity criterion. Distance analyses are based on the overall similarity derived from the values of the parameters describing the objects. The choice of the most adequate distance measure for the data under study is not unique and remains difficult to justify a priori. The way objects are subsequently grouped together (this is called the linkage) is also not uniquely defined. Cladistics uses a specific measure that is based on characters (a trait, a descriptor, an observable, or a property, that can be given at least two states characterizing the evolutionary stages of the object for that character) and compares objects in their evolutionary relationships. Here, the "distance" is an evolutionary cost. Groupings are then made on the basis of shared or inherited characteristics, and are most conveniently represented on an evolutionary tree.

Character-based methods like cladistics are better suited to the study of complex objects in evolution, even though the relative evolutionary costs of the different characters is not easy to assess. Distance-based methods are generally faster and often produce comparable results, but the overall similarity is not always adequate to compare evolving objects. In any case, one has to choose a multivariate method, and the results are generally somewhat different depending on this choice (e.g. Buchanan & Collard (2008)). However, the main goal is to reveal a hidden structure in the data sample, and the relevance of the method is mainly provided by the interpretation and usefulness of the result. In the present paper, the


Julien Menana
email: julien.menana@emn.fr

Sophie Demassey
email: sophie.demassey@emn.fr

Sequencing and Counting with the multicost-regular Constraint

This paper introduces a global constraint encapsulating a regular constraint together with several cumulative costs. It is motivated in the context of personnel scheduling problems, where a schedule meets patterns and occurrence requirements which are intricately bound. The optimization problem underlying the multicost-regular constraint is NP-hard but it admits an efficient Lagrangian relaxation. Hence, we propose a filtering based on this relaxation. The expressiveness and the efficiency of this new constraint is experimented on personnel scheduling benchmark instances with standard work regulations. The comparative empirical results show how multicost-regular can significantly outperform a decomposed model with regular and global-cardinality constraints.

Introduction

Many combinatorial decision problems involve the simultaneous action of sequencing and counting objects, especially in the large class of routing and scheduling problems. In routing, a vehicle visits a sequence of locations following a path in the road network according to some numerical requirements on the whole travelling distance, the time spent, or the vehicle capacity. If only one numerical attribute is specified, finding a route is to solve a shortest/longest path problem. For several attributes, the problem -a Resource Constrained Shortest/Longest Path Problem (RCSPP) -becomes NP-hard. All these numerical requirements may drastically restrict the set of paths in the network which correspond to the actual valid routes. Hence, it is much more efficient to take these requirements into account throughout the search of a path, rather than each separately. Personnel scheduling problems can be treated analogously. Planning a worker schedule is to sequence activities (or shifts) over a time horizon according to many various work regulations, as for example: "a working night is followed by a free morning", "a night shift costs twice as musch as a day shift", "at least 10 days off a month", etc. Hence, a schedule meets both structural requirements -defined as allowed patterns of activities -and numerical requirements -defined as assignment costs or counters -which are intricately bounds. Modelling these requirements individually is itself a hard task, for which the expressiveness and the flexibility of Constraint Programming (CP) is recognized. Modelling these requirements efficiently is still a harder task as it means to aggregate all of them in order to process this set of tied requirements as a whole. By introducing the regular global constraint, Pesant [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF] has proposed an elegant and efficient way to model and to enforce all the pattern requirements together. The allowed patterns are gathered in an acyclic digraph whose paths coincide with the valid sequences of activities. This approach was later extended to optimization constraints soft-regular [START_REF] Van Hoeve | On global warming: Flow-based soft global constraints[END_REF] and cost-regular [START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF] for enforcing bounds on the global cost -a violation cost or any financial cost -of the sequence of assignments. The underlying problem is now to compute shortest and longest paths in the acyclic graph of patterns. The cost-regular constraint was successfully applied to solve real-world personnel scheduling problems under a CP-based column-generation approach [START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF]. Nevertheless, the authors complained about the weak interaction in their CP model between the cost-regular constraint and an external global-cardinality used for modelling occurrence requirements. Actually, with such a decomposition, the support graph of cost-regular maintains many paths which do not satisfy the cardinality constraints. In this paper, we still generalize this approach for handling several cost attributes within one global constraint multicost-regular. Such a constraint allows to reason simultaneously on the sequencing and counting requirements occurring in personnel scheduling problems. As mentioned above, the underlying optimization problem is a RCSPP and it remains NP-hard even when the graph is acyclic. Hence, the filtering algorithm we present achieves a relaxed level of consistency. It is based on the Lagrangian relaxation of the RCSPP following the principle by Sellmann [START_REF] Sellmann | Theoretical foundations of CP-based lagrangian relaxation[END_REF] for Lagrangian relaxation-based filtering. Our implementation of multicost-regular is available in the distribution of the open-source CP solver CHOCO 1 .

The paper is organized as follows. In Section 2, we present the class of regular constraints and provide a theoretical comparison between the pathfinding approach of Pesant [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF] and the decomposition-based approach of Beldiceanu et al. [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF]. We introduce then the new constraint multicost-regular. In Section 3, we introduce the Lagrangian relaxation-based filtering algorithm. In Section 4, we describe a variety of standard work regulations and investigate a systematic way of building one instance of multicost-regular from a set of requirements. In Section 5, comparative empirical results on benchmark instances of personnel scheduling problems are given. They show how multicost-regular can significantly outperform a decomposed model with regular and globalcardinality constraints.

Regular Language Membership Constraints

In this section, we recall the definition of the regular constraint and report on related work, before introducing multicost-regular. First, we recall basic notions of automata theory and introduce notations used throughout this paper:

We consider a non empty set Σ called the alphabet. Elements of Σ are called symbols, sequences of symbols are called words, and sets of words are called languages over Σ. An automaton Π is a directed multigraph (Q, ∆) whose arcs are labelled by the symbols of an alphabet Σ, and where two non-empty subsets of vertices I and A are distinguished. The set Q of vertices is called the set of states of Π, I is the set of initial states, and A is the set of accepting states. The non-empty set ∆ ⊆ Q×Σ ×Q of arcs is called the set of transitions of Π. A word in Σ is said to be accepted by Π if it is the sequence of the arc labels of a path from an initial state to an accepting state in Π. Automaton Π is a deterministic finite automaton (DFA) if ∆ is finite and if it has only one initial state (I = {s}) and no two transitions sharing the same initial extremity and the same label. The language accepted by a FA is a regular language.

Path-Finding and Decomposition: Two Approaches for regular

The regular language membership constraint was introduced by Pesant in [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF]. Given a sequence X = (x 1 , x 2 , ..., x n ) of finite domain variables and a deterministic finite automaton Π = (Q, Σ, ∆, {s}, A), the constraint regular(X, Π) holds iff X is a word of length n over Σ accepted by DFA Π. By definition, the solutions of regular(X, Π) are in one-to-one correspondance with the paths of exactly n arcs connecting s to a vertex in A in the directed multigraph Π. Let δ i ∈ ∆ denote the set of transitions that appears as the i-th arc of such a path, then a value for x i is consistent iff δ i contains a transition labelled by this value.

Coincidently, Pesant [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF] and Beldiceanu et al [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] introduce two orthogonal approaches to achieve GAC on regular (see Figure 1). The approach proposed by Pesant [START_REF] Pesant | A regular language membership constraint for finite sequences of variables[END_REF] is to unfold Π as an acyclic DFA Π n which accepts only the words of length n. By construction, Π n is a layered multigraph with state s in layer 0 (the source), the accepting states A in layer n (the sinks), and where the set of arcs in any layer i coincides with δ i . A breadth-first search allows to maintain the coherence between Π n and the variable domains by pruning the arcs in δ i whose labels are not in the domain of x i , then by pruning the vertices and arcs which are not connected to a source and to a sink. In Beldiceanu et al [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF], a regular is decomposed as n tuple constraints for modelling the sets δ 1 , δ 2 , . . . , δ n . The decomposition introduces state variables q 0 ∈ {s}, q 1 , . . . , q n-1 ∈ Q, q n ∈ A and uses triplet relations defined in extension to enforce GAC on the transition constraints (q i-1 , x i , q i ) ∈ δ i . Such a constraint network being Berge-acyclic, enforcing AC on the decomposition achieves GAC on regular.

In the first approach, a specialized algorithm is defined to maintain all the support paths, while in the second approach, the transitions are modeled with tuple constraints which are directly propagated by the CP solver. The two approaches are orthogonal. Actually, the second model may mimic the specialized algorithm depending on the chosen propagation.

If we assume w.l.o.g. that Σ is the union of the variable domains, then the initial run of Pesant's algorithm for the construction of Π n is performed in O(n|∆|) time and space (with ∆ ≤ |Q||Σ| if Π is a DFA). Incremental filtering is performed with the same worst-case complexity with a forward/backward traversal of Π n . Actually, the complexity of the algorithm relies more on the size |∆ n | of the unfolded automaton Π n rather than on the size |∆| of the specified automaton Π. Note for instance that when the specified automaton Π accepts only words of length n then it is already unfolded (Π = Π n ) and the first run of the algorithm is in O(|∆|). In practice, as in our experiments (Section 5), Π n can even be much smaller than Π, meaning that many accepting states in Π cannot be reached in exactly n transitions. The incremental filtering is performed in

O(|∆ n |) time with, in such a case, |∆ n | ≪ n|∆|.
regular is a very expressive constraint. It is useful to model pattern constraints arising in many planning problems, but also to reformulate other global constraints [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] or to model tuples defined in extension. An other application of regular is to model a sliding constraint: recently, Bessière et al. [START_REF] Bessière | Reformulating global constraints: The SLIDE and REGULAR constraints[END_REF] have introduced the slide meta constraint. In its more general form, slide takes as arguments a matrix of variables Y of size n × p and a constraint C of arity pk with k ≤ n. slide(Y, C) holds if and only of C(y 1 i+1 , . . . , y p i+1 , . . . , y 1 i+k , . . . , y p i+k ) holds for 0 ≤ i ≤ n-k. Using the decomposition proposed in [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF], regular(X, Π) can be reformulated as slide([Q, X], C ∆ ), where Q is the sequence of state variables and C ∆ is the transition constraint C ∆ (q, x, q ′ , x ′ ) ≡ (q, x, q ′ ) ∈ ∆. Conversely [START_REF] Bessière | Reformulating global constraints: The SLIDE and REGULAR constraints[END_REF], a slide constraint can be reformulated as a regular but it may require to enumerate all valid tuples for C. This reformulation can however be useful in the context of planning (especially for car sequencing) to model a sliding cardinality constraint also known as sequence. Even if powerful specialized algorithms exist for this constraint (see e.g. [START_REF] Maher | Flow-Based Propagators for the sequence and Related Global Constraints[END_REF]), the automaton resulting from the reformulation can be integrated with other pattern requirements as we will show in Section 4. Finally, one should notice the work (see e.g. [START_REF] Kadioglu | Efficient context-free grammar constraints[END_REF]) related to contextfree grammar constraints. Though, most of the rules encountered in personnel scheduling can be described using regular languages.

Maintaining Patterns with Cumulative Costs and Cardinalities

Personnel scheduling problems are usually defined as optimization problems. Most often, the criterion to optimize is a cumulative cost, i.e. the sum of costs associated to each assignment of a worker to a given activity at a given time.

Such a cost has several meanings: it can model a financial cost, a preference, or a value occurrence. Now, designing a valid schedule for one worker is to enforce the sequence of assignments to comply with a given pattern while ensuring that the total cost of the assignments is bounded. This can be specified by means of a cost-regular constraint [START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF]. Given c = (c ia ) i∈ [1..n]×a∈Σ a matrix of real assignment costs and z ∈ [z, z] a bounded variable (z, z ∈ R), cost-regular(X, z, Π, c) holds iff regular(X, Π) holds and n i=1 c ixi = z. Note that it has the knapsack constraint [START_REF] Trick | A dynamic programming approach for consistency and propagation for knapsack constraints[END_REF] as a special case and that, unless P = N P , one can enforce GAC on a knapsack constraint at best in pseudo-polynomial time, i.e. the run time is polynomial in the values of the bounds of z. As a consequence, enforcing GAC on cost-regular is NP-hard.

The definition of cost-regular reveals a natural decomposition as a regular constraint channeled to a knapsack constraint. Actually, it is equivalent to the decomposition proposed by Beldiceanu et al. [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] when dealing with one cumulative2 cost: cost variables k i are now associated to the previous state variables q i , with k 0 = 0 and k n = z, and several arithmetic and element constraints model the knapsack and channeling constraints. In short, this formulation can be rewritten as slide

([Q, X, K], C c ∆ ), with C c ∆ (q i-1 , x i-1 , k i-1 , q i , x i , k i ) ≡ (q i-1 , x i-1 , q i ) ∈ ∆∧k i = k i-1 +c ixi .
Depending on the size of the domains of the cost variables, GAC can be enforced on knapsack in reasonable time. However, even in this case, since the constraint hypergraph of the decomposed model is no longer Berge-acyclic but α-acyclic, one has to enforce pairwise-consistency on the shared variables -a pair (q i , k i ) of state and cost variables -of the transition constraints in order to achieve GAC. A similar option proposed for slide [START_REF] Bessière | Reformulating global constraints: The SLIDE and REGULAR constraints[END_REF] is to enforce AC on the dual encoding of the hypergraph of the C c ∆ constraints, but again it requires to explicit all the support tuples and then, it may be of no practical use.

The filtering algorithm presented in [START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF] for cost-regular is a slight adaptation3 of Pesant's algorithm for regular. It is based on the computation of shortest and longest paths in the unfolded graph Π n valued by the transition costs. To each vertex (i, q) in any layer i of Π n are associated two bounded cost variables k - iq and k + iq modelling the lengths of the paths respectively from layer 0 to (i, q) and from (i, q) to layer n. The cost variables can trivially be initialized during the construction of Π n : k - iq in the forward phase and k + iq in the backward phase. The bounds of variable z are then pruned according to the condition z ⊆ k + 0s . Conversely, an arc ((i -1, q), a, (i, q ′ )) ∈ δ i can be removed whenever:

k - (i-1)q + c ia + k + iq ′ > z or k - (i-1)q + c ia + k + iq ′ > z.
As graph Π n is acyclic, maintaining the cost variables, i.e. shortest and longest paths, can be performed by breadth-first traversal with the same time complexity O(|∆ n |) than for maintaining the connexity of the graph in regular.

As said before, this algorithm achieves a hybrid level of consistency on cost-regular. As a matter of fact, it enforces a sort of pairwise-consistency on the decomposed model between each state variable and the bounds of the associated cost variable, according to the relation q i = (i, q) ⇐⇒ k i = k - iq . Hence, it dominates the decomposed model knapsack∧regular when only Bound Consistency is enforced on the cost variables. Otherwise, if AC is enforced on knapsack then the two approaches are incomparable as show the two examples depicted in Figures 2 and3.

3 s 1 2 a [1] b [0] a [1] [0,1] b [0] x1 ∈ {// a, b}, x2 ∈ {// a, b}, z ∈ [0, 1]. (q0, x1, q1, j1) ∈ {(s, a, 1, 1), (s, b, 2, 2)}, element(k1, j1, (k0 + 1, k0)), (q1, x2, q2, j2) ∈ {(1, a, 3, 1), (2, b, 3, 2)}, element(z, j2, (k1 + 1, k1)), q0 ∈ {s}, q1 ∈ {1, 2}, q2 ∈ {3}, k0 ∈ {0}, k1 ∈ {0, 1}, x1 ∈ {a, b}, x2 ∈ {a, b}, z ∈ {0, 1}.
[2,2] a [0] b [2] b [2] c [1] a [0] 2 1 x1 ∈ {a, b, c}, x2 ∈ {a, b}, z ∈ [2, 2].
(x1, j1) ∈ {(a, 1),/////// (c, 2), (b, 3)}, element(k1, j1, (k0,///////// k0 + 1, k0 + 2)), (x2, j2) ∈ {(a, 1), (b, 2)}, element(z, j2, (k1, k1 + 2)), k0 ∈ {0}, k1 ∈ {0, // 1, 2}, x1 ∈ {a, b,/ c}, x2 ∈ {a, b}, z ∈ {2}.

The multicost-regular Constraint

A natural generalization of cost-regular is to handle several cumulative costs: given a vector Z = (z 0 , ..., z R ) of bounded variables and c = (c r ia )

r∈[0..R] i∈[1.
.n],a∈Σ a matrix of assignment costs, multicost-regular(X, Z, Π, c) holds if and only if regular(X, Π) holds and n i=1 c r ixi = z r for all 0 ≤ r ≤ R. Such a generalization has an important motivation in the context of personnel scheduling. Actually, apart a financial cost and pattern restrictions, an individual schedule is usually subject to a global-cardinality constraint bounding the number of occurrences of each value in the sequence. These bounds can drastically restrict the language on which the schedule is defined. Hence, it could be convenient to tackle them within the regular constraint in order to reduce the support graph. As a generalization of cost-regular or of the global-sequencing constraint [START_REF] Régin | A filtering algorithm for global sequencing constraints[END_REF], we cannot hope to achieve GAC in polynomial time here. Note that the model by Beldiceanu et al [START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] -and similarly the slide constraint -was also proposed for dealing with several costs but again, it amounts to decompose as a regular constraint channeled with one knapsack constraint for each cost.

Hence, we ought to exploit the structure of the support graph of Π n to get a good relaxed propagation for multicost-regular. The optimization problems underlying cost-regular were shortest and longest path problems in Π n . The optimization problems underlying multicost-regular are now the Resource Constrained Shortest and Longest Path Problems (RCSPP and RCLPP) in Π n . The RCSPP (resp. RCLPP) is to find the shortest (resp. longest) path between a source and a sink in a valued directed graph, such that the quantities of resources accumulated on the arcs do not exceed some limits. Even with one resource on acyclic digraphs, this problem is known to be NP-hard [START_REF] Handler | A dual algorithm for the restricted shortest path problem[END_REF]. Two approaches are most often used to solve RCSPP [START_REF] Handler | A dual algorithm for the restricted shortest path problem[END_REF]: dynamic programming and Lagrangian relaxation. Dynamic programming-based methods extend the usual short path algorithms by recording the costs over every dimension at each node of the graph. As in cost-regular, this could easily be adapted for filtering by converting these cost labels as cost variables but it would make the algorithm memory expensive. Instead, we investigate a Lagrangian relaxation approach, which can also easily be adapted for filtering from the cost-regular algorithm without memory overhead.

A Lagrangian Relaxation-Based Filtering Algorithm

Sellmann [START_REF] Sellmann | Theoretical foundations of CP-based lagrangian relaxation[END_REF] laid the foundation for using the Lagrangian relaxation of a linear program to provide a cost-based filtering for a minimization or maximization constraint. We apply this principle to the RCSPP/RCLPP for filtering multicost-regular. The resulting algorithm is a simple iterative scheme where filtering is performed by cost-regular on Π n for different aggregated cost functions. In this section, we present the usual Lagrangian relaxation model for the RCSPP and explain how to solve it using a subgradient algorithm. Then, we show how to adapt it for filtering multicost-regular.

Lagrangian Relaxation for the RCSPP. Consider a directed graph G = (V, E, c) with source s and sink t, and resources (R 1 , ..., R R ). For each resource 1 ≤ r ≤ R, let z r (resp. z r ) denote the maximum (resp. minimum 4 ) capacity of a path over the resource r, and c r ij denote the consumption of resource r on arc (i, j) ∈ E. A binary linear programming formulation for the RCSPP is as follows:

min (i,j)∈E c ij x ij (1)
s.t. z r ≤ (i,j)∈E c r ij x ij ≤ z r ∀r ∈ [1..R] (2) j∈V x ij - j∈V x ji =    1 if i = s, -1 if i = t, 0 otherwise. ∀i ∈ V (3)
x ij ∈ {0, 1} ∀(i, j) ∈ E. (4)
In this model, a binary decision variable x ij defines whether arc (i, j) belongs to a solution path. Constraints (2) are the resource constraints and Constraints (3) are the usual path constraints.

Lagrangian relaxation consists in dropping "complicating constraints" and adding them to the objective function with a violation penalty cost u ≥ 0, called the Lagrangian multipliers. The resulting program is called the Lagrangian subproblem with parameter u and it is a relaxation of the original problem. Solving the Lagrangian dual is to find the multipliers u ≥ 0 which gives the best relaxation, i.e. the maximal lower bound.

The complicating constraints of the RCSPP are the 2R resource constraints (2). Indeed, relaxing these constraints leads to a shortest path problem, that can be solved in polynomial time. Let P denote the set of solutions x ∈ {0, 1} E satisfying Constraints (3). P defines the set of paths from s to t in G. The Lagrangian subproblem with given multipliers u = (u -, u + ) ∈ R 2R + is:

SP (u) : f (u) = min x∈P cx + R r=1 u r + (c r x -z r ) - R r=1 u r -(c r x -z r ) (5)
An optimal solution x u for SP (u) is then a shortest path in graph G(u) = (V, E, c(u)) where:

c(u) = c + R r=1 (u r + -u r -)c r , κ u = n r=1 (u r -z r -u r + z r ) and f (u) = c(u)x u + κ u . ( 6
)
Solving the Lagrangian Dual. The Lagrangian dual problem is to find the best lower bound f (u), i.e. to maximize the piecewise linear concave function f :

LD : f LD = max u∈R 2R + f (u) (7)
Several algorithms exist to solve the Lagrangian dual. In our approach, we consider the subgradient algorithm [START_REF] Shor | Minimization methods for non-differentiable functions[END_REF] as it is rather easy to implement and it does not require the use of a linear solver. The subgradient algorithm iteratively solves one subproblem SP (u) for different values of u. Starting from an arbitrary value, the position u is updated at each iteration by moving in the direction of a supergradient Γ of f with a given step length µ: u p+1 = max{u p + µ p Γ (u p ), 0}. There exist many ways to choose the step lengths for guaranteeing the convergence of the subgradient algorithm towards f LD (see e.g. [START_REF] Boyd | Subgradient methods[END_REF]). In our implementation, we use a standard step length µ p = µ 0 ǫ p with µ 0 and ǫ < 1 "sufficiently" large (we have empirically fixed µ 0 = 10 and ǫ = 0.8). For the supergradient, solving SP (u) returns an optimum solution x u ∈ P and Γ (u) is computed as:

Γ (u) = ((c r x u -z r ) r∈[1..R] , (z r -c r x u ) r∈[1..R] ).
From Lagrangian Relaxation to Filtering. The key idea of Lagrangian relaxation-based filtering, as stated in [START_REF] Sellmann | Theoretical foundations of CP-based lagrangian relaxation[END_REF], is that if a value is proved to be inconsistent in at least one Lagrangian subproblem then it is inconsistent in the original problem:

Theorem 1. (i) Let P be a minimization linear program with optimum value f * ≤ +∞, z ≤ +∞ be an upper bound for P , and SP (u) be any Lagrangian subproblem of P , with optimum value

f (u) * ≤ +∞. If f (u) > z then f * > z.
(ii) Let x be a variable of P and v a value in its domain. Consider P x=v (resp. SP (u) x=v ) the restriction of P (resp. SP (u)) to the set of solutions satisfy-

ing x = v and let f * x=v ≤ +∞ (resp. f (u) x=v ≤ +∞) its optimum value. If f (u) x=v > z then f * x=v > z.
Proof. Statement (i) of Theorem 1 is straightforward, since SP (u) is a relaxation for P , then f (u) ≤ f * . Statement (ii) arises from (i) and from the fact that, adding a constraint x = v within P and applying Lagrangian relaxation, or applying Lagrangian relaxation and then adding constraint x = v to each subproblem, result in the same formulation.

The mapping between multicost-regular(X, Z, Π, c), with |Z| = R+1 and an instance of the RCSPP (resp. RCLPP) is as follows: We single out one cost variable, for instance z 0 , and create R resources, one for each other cost variable. The graph G = (Π n , c 0 ) is considered. A feasible solution of the RCSPP (resp. RCLPP) is a path in Π n from the source (in layer 0) to a sink (in layer n) that consumes on each resource 1 ≤ r ≤ R is at least z r and at most z r . Furthermore, we want to enforce an upper bound z 0 on the minimal value for the RCSPP (resp. a lower bound z 0 on the maximal value for the RCLPP). The arcs of Π n are in one-to-one correspondance with the binary variables in the linear model of these two instances.

Consider a Lagrangian subproblem SP (u) of the RCSPP instance (the approach is symmetric for the maximization instance of RCLPP) . We show that a slight modification of the cost-regular algorithm allows to solve SP (u) but also to prune arcs of Π n according to Theorem 1 and to shrink the lower bound z 0 . The algorithm starts by updating the costs on the graph Π n with c 0 (u), as defined in [START_REF] Bessière | Reformulating global constraints: The SLIDE and REGULAR constraints[END_REF] and then by computing, at each node (i, q), the shortest path k - iq from layer 0 and the shortest path k + iq to layer n. We get the optimum value

f (u) = k + 0s + κ u .
As it is a lower bound for z 0 , one can eventually update this lower bound as z 0 = max{f (u), z 0 }. Then, by a traversal of Π n , we remove each arc ((i -1, q), a, (i, q ′ )) ∈ δ i such that k - (i-1)q + c 0 (u) ia + k + iq ′ > z 0κ u . The global filtering algorithm we developed for multicost-regular is as follows: starting from u = 0, a subgradient algorithm guides the choice of the Lagrangian subproblems to which the above cost-filtering algorithm is applied. The number of iterations for the subgradient algorithm is limited to 20 (it usually terminates far before). The subgradient algorithm is first applied to the minimization problem (RCSPP) then to the maximization problem (RCLPP). As a final step, we run the original cost-regular algorithm on each of the cost variables to shrink their bounds (by the way, it could deduce new arcs to filter, but it did not happen in our experiments). Note that due to the parameter dependancy of the subgradient algorithm, the propagation algorithm is not monotonic.

Modelling Personnel Scheduling Problems

In this section we show how to model standard work regulations arising in Personnel Scheduling Problems (PSP) as one instance of the multicost-regular constraint. The purpose is to emphasize the ease of modelling with such a constraint and also to derive a systematic way of modelling PSP.

Standard work regulations

In PSP, many kinds of work regulations can be encountered, however, we can categorize most of them as rules enforcing either regular patterns, fixed cardinalities or sliding cardinalities.

To illustrate those categories, we consider a 7 days schedule and 3 activities: night shift (N), day shift (D) or rest shift (R). For example: R R D D N R D Regular patterns can be modelled directly as a DFA. For instance the rule "a night shift is followed by a rest" is depicted in Figure 4 (A). The rules can either be given as forbidden patterns or allowed patterns. In the first case, one just need to build the complement automaton. (A) a regular pattern example. Fixed cardinality rules bound the number of occurrences of an activity or a set of activities over a fixed subsequence of time slots. Such a rule can be modelled within an automaton or using counters. For example, the rule "at least 1 and at most 3 day shifts each week" can easily be modelled as the DFA depicted in Figure 4 (B). Taking a look at this automaton, we can see the initial state has been split into 3 different states that represent the maximum number of D transitions that can be taken. Such a formulation can be an issue when the maximum occurrence number increases. In this case, using a counter is more suitable, as we only need to create a new cost variable z r ∈ [1, 3] with c r ij = 1 ⇐⇒ 1 ≤ i ≤ 7 and j = D. More generally, one can also encounter cardinality rules over patterns. This also can be managed by means of a cost. One has to isolate the pattern within the automaton describing all the feasible schedules, then to price transitions entering it to 1.

Sliding constraints can be modelled as a DFA using the reformulation stated in [START_REF] Bessière | Reformulating global constraints: The SLIDE and REGULAR constraints[END_REF]. However, the width of the sliding sequence should not be too large as the reformulation requires to explicit all the feasible tuples of the constraint to slide. This is often the case in PSP or also in car sequencing problems.

Systematic multicost-regular Generation

A formalism to describe Personnel Scheduling Problems has been proposed in [14]. The set of predefined XML markups allows to specify a large scope of PSP. In order to automatically generate a CP model based on multicost-regular from such specifications, we developed a framework capable of interpreting those XML files. In a first step, we bounded each markup associated to a work regulation to one of the 3 categories described above. Hence, for each rule of a given PSP instance, we automatically generate either an automaton or a counter depending on the rule category. For instance, the forbidden pattern "no day shift just after a night shift" is defined in the xml file as <Pattern weight="1350"><Shift>N</Shift><Shift>D</Shift></Pattern> and is automatically turned into its equivalent regular expression (D|N |R) * N D(D|N |R) * . We use a java library for automata5 in order to create a DFA from a regular expression and to operate on the set of generated DFA. We use the opposite, the intersection and the minimization operations to build an unique DFA. Once the DFA is built, we treat the rules that engender counters, and generate a multicost-regular instance for each employee. Last, we treat the transversal constraints and include them in the CP model. For example, cover requirements are turned into global-cardinality constraints. Note that we were not able to deal with two kinds of specifications: some rule violation penalties that the multicost-regular cannot model and the pattern cardinality rules that we do not yet know how to automatize the reformulation.

Two Personnel Scheduling Cases

We first tackled the GPost [14] problem. This PSP consists in building a valid schedule of 28 days for eight employees. Each day, an employee has to be assigned to a Day, Night, or Rest Shift. Each employee is bound to a (Fulltime or Parttime) contract defining regular pattern and cardinality rules. Regular pattern rules are: "free days period should last at least two days", "consecutive working week-ends are limited" and "given shift sequences are not allowed". Using the automatic modelling method we presented earlier, we build a DFA for each kind of contract. Cardinality rules are: "a maximum number of worked days in the 28 days period is to be worked", "the amount of certain shifts in a schedule is limited" and "the number of working days per week is bounded". Cover requirements and employee availabilities are also modelled. The softness specification on rules has been ignored as well as the first pattern rule to avoid infeasibility.

The second case study is based on the generated benchmark set brought by Demassey et al. [START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF]. The work regulations arise from a real-world personnel scheduling problem. The goal is to build only one schedule for a day consisting of 96 fifteen minutes time slots. Each slot is assigned either a working activity, a break, a lunch or a rest. Each possible assignment carries a given cost. The purpose is to find a schedule of minimum cost meeting all the work regulations. As for the previous PSP, we can identify regular pattern work regulations: "A working activity lasts at least 1 hour", "Different work activities are separated by a break or a lunch", "Break, lunch and rest shifts cannot be consecutives", "Rest shifts are at the beginning or at the end of the day", and "A break lasts 15 minutes". And also fixed cardinality regulations with: "At least 1 and at most 2 breaks a day", "At most one lunch a day" and "Between 3 and 8 hours of work activities a day". In addition to those work regulations, some activities are not allowed to be performed during some period. These rules are trivial to model with unary constraints.

Experiments

Experimentations were run on an Intel Core 2 Duo 2Ghz processor with 2048MB of RAM running OS X. The two PSP problems were solved using the Java constraint library CHOCO with default value selection heuristics -min value.

On the Size of the Automaton

As explained in Section 2.1, the filtering algorithm complexity of the regular constraints depends on the size of the specified automaton. Thus it would seem natural that processing a big automaton is not a good idea. However, practical results points out two important facts. First of all, the operations we run for automatically building a DFA from several rules tends to generate partially unfolded DFA (by intersection) and to reduce the number of redundant states which lie in the same layers (by minimization). Hence, the unfolded automaton generated during the forward phase at the initialization of the constraint can even be smaller than the specified automaton Π. Secondly, pruning during the backward phase may produce an even smaller automaton Π n as many accepting states cannot be reached in a given number n of transitions. Table 1 shows the number of nodes and arcs of the different automata during the construction of the multicost-regular constraints for the GPost problem: the sum of the DFAs generated for each rule, the DFA Π after intersection and minimization, the unfolded DFA after the forward and backward phases.

Comparative Experiments

The previous section showed the ease of modelling with multicost-regular. However, there would be no point in defining such a constraint if the solving was badly impacted. We then conduct experiments for comparing our algorithm with a decomposed model consisting of a regular (or cost-regular for optimization) channeled to a global-cardinality constraint (gcc).

Table 2 presents the computational results on the GPost instance. The models include 8 multicost-regular or 8 regular and gcc (for each employee) bound together by 28 transversal gcc (for each day). In the Table, the first row corresponds to the problem without the sliding rule over the maximum number of consecutive working week-ends. In the second row, this constraint was included. We tried various variable selection heuristics but found out assigning variables along the days gave the best results as it allows the constraint solver to deal with the transversal gcc more efficiently. Both models lead to the same solutions. Actually, the average time spent on each node is much bigger using multicost-regular. However, due to better filtering capabilities, the size of the search tree and the runtime to find a feasible solution are significantly decreased.

Our second experiment tested the scalability of multicost-regular (MCR) against cost-regular∧gcc (CR) on the optimization problem defined in Section 4.3. The models do not contain any other constraint. However the decomposed model CR requires additional channeling variables. Table 3 presents the results on a benchmark set made of 110 instances. The number n of working activities varies between 1 and 50. The assignment costs were randomly generated. We tested different variable selection heuristics and kept the best one for each model. Note that the results of the CR model are more impacted by the heuristic.

The first columns in the table show that with the MCR model, we were able to solve all instances (Column #) in less than 15 seconds for the biggest ones (Column t). The average number of backtracks (Column bt) remains stable and low as n increases. On the contrary the CR model is impacted a lot as shown in the next columns. Indeed, as the initial underlying graph becomes bigger it contains more and more paths violating the cardinality constraints. Those paths are not discarded by cost-regular. Some instances with more than 8 activities could not be solved within the given 30 minutes (Column #). Considering only solved instances, the running time (t) and the number of backtracks (bt) are always much higher than the MCR model based results. Regarding unsolved instances, the best found solution within 30 minutes is rarely optimal (Column # opt), and the average gap (Column ∆) is up to 6% for 40 activities.

Conclusion

In this paper, we introduce the multicost-regular global constraint and provide a simple implementation of Lagrangian relaxation-based filtering for it. Experimentations on benchmark instances of personnel scheduling problems show the efficiency and the scalability of this constraint compared to a decomposed model dealing with pattern requirements and cardinality requirements separately. Furthermore, we investigate a systematic way to build an instance of multicost-regular from a given set of standard work regulations. In future works, we ought to get a fully systematic system linked to the CHOCO solver for modelling and solving a larger variety of personnel scheduling and rostering problems.

Fig. 1 .

 1 Fig. 1. Consider the DFA depicted above applied to X ∈ {a, b}×{a}×{a, b}×{b}. The unfolded automaton of regular is depicted on the left and the decomposed model on the right. The dashed transitions are discarded in both models.

Fig. 2 .

 2 Fig. 2. Consider the depicted DFA with costs in brackets applied to X = (x1, x2) ∈ {a, b} × {a, b} and z ∈ [0, 1]. The cost-regular algorithm (on the left) discards the dashed transitions and hence achieves GAC. The decomposed model (on the right) is arc-consistent but not globally consistent.

  s

Fig. 3 .

 3 Fig. 3. Consider now the depicted DFA applied to X = (x1, x2) ∈ {a, b, c} × {a, b} and z ∈ [2, 2]. Enforcing AC on the decomposed model (on the right) achieves GAC. The cost-regular algorithm (on the left) does not achieve GAC since the minimum and maximum paths traversing arc x1 = c are consistent with the bounds on z.

Fig. 4 .

 4 Fig. 4. Examples of automata representing work regulations.

Table 1 .

 1 Illustration of graph reduction during presolving.

	Contract Count	sum of DFAs Π Forward Backward (Πn)
	Fulltime	# Nodes # Arcs	5782 682 40402 4768	411 1191	230 400
	Parttime	# Nodes # Arcs	4401 385 30729 2689	791 2280	421 681

Table 2 .

 2 GPost problem results

		multicost-regular regular ∧ gcc
	WE regulation Time (s)	# Fails Time (s) # Fails
	no	1.94	24	12.6 68035
	yes	16.0	1576	449.2 2867381

Table 3 .

 3 Shift generation results

		MCR model				CR model
			Solved			Solved			Time out
	n	#	t	bt	#	t	bt	#	∆	# opt
	1	10	0.6	49	10	1.1	292	0	-	-
	2	10	0.8	54	10	2.4	539	0	-	-
	4	10	1.5	65	10	13.5	1638	0	-	-
	6	10	1.6	44	10	53.6	4283	0	-	-
	8	10	2.1	51	9	209.2	5132	1	3.5%	0
	10	10	2.4	58	7	283.5	6965	3	4.6%	0
	15	10	3.8	59	6	283.9	4026	4	4.7%	1
	20	10	4.9	49	6	311.8	4135	4	4.2%	1
	30	10	6.9	51	1	313.0	4303	9	3.1%	1
	40	10	13.4	68	0	-	-	10	6.1%	0
	50	10	14.4	51	1	486.0	1406	9	5.0%	1

http://choco.emn.fr/

The model in[START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] can deal not only with sum but also with various arithmetic functions on costs, but no example of use is provided.

[START_REF] Demassey | A cost-regular based hybrid column generation approach[END_REF] Previously, the algorithm was partially -for minimization only -applied to the special case soft-regular[hamming] in[START_REF] Beldiceanu | Reformulation of Global Constraints Based on Constraint Checkers[END_REF] and in[START_REF] Van Hoeve | On global warming: Flow-based soft global constraints[END_REF].

in the original definition of RCSPP, there is no lower bound on the capacity: zr

http://www.brics.dk/automaton/

Acknowledgements

We thank Mats Carlsson for pointing out the example illustrating the propagation issue of the cost-regular algorithm. We also thank Christian Schulte for its insightful comments on the paper and its numerous ideas to improve the constraint.

Data Sets and Benchmarks: http://www.cs.nott.ac.uk/~tec/NRP/