Keywords: Optimal shipment size, freight transportation economics, inventory theory, empirical analysis, EOQ model

Deciding shipment size is important in freight transport: it depends on the logistical imperatives of shippers and the technical possibilities of carriers. Shipment size choice is also closely related to transportation mode; it is therefore important from a public policy perspective.

The theory of optimal shipment size and mode choice is robust. There are many inventory-theoretical models of optimal shipment size, applied by shippers in operational contexts. However, none of them has been validated empirically over a large and heterogeneous population of shipments, so that they are virtually useless for modeling freight transportation demand. This is due in particular to the lack of adequate data.

In this paper, the simple Economic Order Quantity (EOQ) model is assessed empirically on a national scale, over a heterogeneous population of shipments. Using the French ECHO database, which notably observes commodity flow rates between shippers and receivers, the EOQ shipment size specification is estimated. The validity of the EOQ model is confirmed. In addition, the dominant role of the commodity flow rate between the shipper and the receiver, and of commodity value density is revealed. The relationship between mode choice and shipment size is also highlighted.

INTRODUCTION

Freight transportation demand models are largely inspired by passenger transportation demand models [START_REF] De Jong | National and international freight transport models: an overview and ideas for future development[END_REF]. Their architecture is generally based closely on the classic four stages of passenger transportation demand modeling [START_REF] Ortuzar | Modelling transport[END_REF] i.e. generation, distribution, mode choice, and assignment. In particular, classical passenger mode choice models are generally discrete choice models, where the value of a mode alternative to a passenger is represented by additive random utility functions (3; 4; 5). These utility functions consist of variables associated to the individuals and of variables associated to the transport alternatives (typically their rates, travel time, and reliability). Freight mode choice models generally follow the same principles.

These models are relatively efficient empirically, but they have shortcomings. Probably one of the most significant is the seemingly total neglect of the role of logistics. This is an issue for at least three reasons. First because it suggests a theoretical flaw: for example, it is impossible to establish a formal relationship between the preferences of a shipper regarding freight transport and its own logistical environment (customer preferences, position in the supply chain, etc.) Second, many transport policy instruments fall outside the scope of these models (warehouse location policies, taxes on warehousing, etc.) Third, an important variable lacks: shipment size.

Shipment size plays an important role in freight transportation. For shippers, sending ten shipments of a thousand tons a year is entirely different from sending ten thousands oneton shipments. For carriers, different shipment sizes mean different vehicles and organizations. This applies especially to the choice of transportation mode.

Theoretical models of optimal shipment size and mode choice have been around for a long time. They are most often based on the century-old Economic Order Quantity model. Initially designed to optimize production in certain contexts, it applies easily to the context of freight transport. However, its validity for a large and heterogeneous population of shippers has never yet been assessed. This is mainly due to the lack of an adequate database. Indeed, an important explanatory variable in the EOQ model is the commodity flow rate between the shipper and the receiver. This variable is not usually recorded. The French shipment database ECHO is one exception. This paper uses that database in order to reach an empirical assessment of the Economic Order Quantity model for shipment size choice in freight transport.

The discussion proceeds as follows. First, Section 2 presents a review of the literature mainly focused on the choice of shipment size in freight transport. Section 3 recalls the theoretical underpinning and formulas of the EOQ model. Section 4 then presents the ECHO database and the estimation methodology. Section 5 discusses the results of that estimation. Section 6 presents the estimation of an extended EOQ model, which gives additional insight into the structure of transportation costs and their impact on the choice of shipment size. The paper is concluded in Section 6.

LITERATURE REVIEW

Freight mode choice models are generally based on microeconomic, often disaggregate models of the interaction of shippers and carriers. Initial investigations on freight mode choice were focused on direct transport costs (6; 7). Later studies concluded that indirect costs played an important role (8; 9; and also 10). Some models were based on total cost minimization and linear programming [START_REF] Casavant | Multimodal transportation and impacts of policy: grain transportation model[END_REF]. With the parallel development of discrete choice modeling techniques [START_REF] Ben-Akiva | Discrete Choice Analysis: Theory and Application to Travel Demand[END_REF], this led to a now dominant model architecture where each modal alternative is represented by a utility (or generalized cost) and the choice between alternatives by a probabilistic model, such as the multinomial logit and its refinements (13, 14; for a comprehensive review of recent freight models, refer to 15; for a review of some recent discrete choice model specifications in freight transport, refer to 16). All these models have in common that shipment size is absent.

The importance of shipment size was identified long ago. It was soon found that shipment size is an empirically important explanatory variable in mode choice (17; 18). [START_REF] Winston | A disaggregate model of the demand for intercity freight transportation[END_REF] designed and estimated a behavioral model where the transportation mode is assumed to be chosen by shippers (or receivers) for given, exogenous shipment sizes. Another example is [START_REF] Jiang | Freight demand characteristics and mode choice: an analysis of the results of modeling with disaggregate revealed preference data[END_REF], where the authors estimated a disaggregate model of mode choice using shipment characteristics as explanatory variables, on the basis of the 1988 French shipment database Enquête Chargeurs.

However, it was soon understood that shipment size depends as much on mode choice as mode choice depends on shipment size. This was confirmed empirically using joint estimation techniques [START_REF] Abdelwahab | Modeling the demand for freight transport: new approach[END_REF]. In this study, a switching regression model was estimated. The shipment sizes conditional to each mode were specified as linear combinations of exogenous parameters. Then, a latent variable, specified as a linear combination of these parameters and of the two conditional shipment sizes, was used itself to predict the transportation mode, on the basis of its sign. A similar approach, yielding similar qualitative results, is presented in [START_REF] Holguin-Veras | Revealed preference analysis of commercial vehicle choice process[END_REF]. In [START_REF] Windisch | A disaggregate freight transport model of transport chain and shipment size choice[END_REF], shipment size was modeled as a discrete variable, so that the discrete choice modeling toolbox could be applied straightforwardly. The interdependency of mode and shipment size was confirmed again. All those approaches led to the same empirical conclusion. Nevertheless, they gave limited insight on its underlying microeconomic drivers.

Theoretical models of optimal shipment size are more than a century old. They are part of inventory theory. The first one is the EOQ -Economic Order Quantity -model (24; 25). It was initially derived to provide the optimal batch size for a production chain, but it applies straightforwardly to freight transport. According to this model (described in more detail in the next section) the optimal shipment size is the result of a trade-off between transportation costs and inventory costs. It can be extended to take mode choice into account as well, and even other variables. This is the approach chosen in the very important paper [START_REF] Baumol | An inventory theoretic model of freight transport demand[END_REF], which derived a theoretical, disaggregate model of a combined decision on transportation mode, shipment size, and safety inventory at destination. This paper was followed by a series of contributions, of which we will mention a small, significant sample. The influence of the shape of freight rates on choice of shipment size and transportation mode was investigated in [START_REF] Langley | The inclusion of transportation costs in inventory models: some considerations[END_REF]. The model introduced in (26) was simplified to represent the choice of shipment size and type of transportation operation by road, between parcel, less-than-truckload and truckload shipments [START_REF] Hall | Dependence between shipment size and mode in freight transportation[END_REF]. A dynamic inventory-theoretical model is derived and estimated in [START_REF] Mcfadden | Joint estimation of freight transportation decisions under nonrandom sampling[END_REF]. However, the econometric specification was only loosely based on the theoretical developments from which it purportedly derives. In [START_REF] Mccann | A proof of the relationship between optimal vehicle size, haulage length and the structure of distance-transport costs[END_REF], the EOQ model formed the basis of a theoretical model of optimal ship size. A significant current goal of research in freight transportation modeling is to build a fully-fledged freight transportation demand model with an explicit representation of shipment size. Such an attempt is described in [START_REF] De Jong | A micro-simulation model of shipment size and transport chain choice[END_REF]. Finally, the overall status of inventory theory in freight transportation modeling is discussed in (32; 33; 34).

Despite their theoretical appeal, these models all lack solid empirical validation. This is mainly due to the absence of adequate data. For example, to estimate the EOQ model, one needs to observe shipment sizes as well as shipper-receiver commodity flow rates. While there are databases where shipment size is available -such as the US or Swedish Commodity Flow Surveys (35; 36), the commodity flow rate between shipper and receiver is overwhelmingly ignored. Some authors have been able to access relevant datasets, but these are restricted to very small samples (no more than two firms). In such cases, estimating an inventory-theoretical model may determine whether a particular logistical doctrine holds for those firms, but this does not constitute evidence that it is valid at an aggregate scale (4; 37; 38). Experimental economics can partially address this issue, but not as satisfyingly as real world data [START_REF] Holguin-Veras | An experimental economics investigation of shipper-carrier interactions in the choice of mode and shipment size in freight transport[END_REF]. Attempts to replace the shipper-receiver commodity flow with proxies in a large-scale spatialized model have encountered difficulties [START_REF] De Jong | A micro-simulation model of shipment size and transport chain choice[END_REF]. In consequence, no functional spatialized freight transportation model currently includes shipment size as an endogenous variable [START_REF] Chow | State-of-the-art of freight forecast modelling: lessons learned and the road ahead[END_REF].

An adequate disaggregate shipment database is required to assess the empirical validity of the EOQ model. In particular, it should describe not only the characteristics of shipments, and of the way they are transported, but also the logistical contexts in which these transportation operations take place. One database provides such information: the French shipper survey ECHO (Enquête CHargeurs et Opérateurs de transport, shippers and transport operators survey), (41; 42). This is the database used in this paper.

THE ECONOMIC ORDER QUANTITY MODEL

The recent discussions on freight transportation modeling insist that supply and demand should be clearly distinguished (43; 44). In this paper, we follow [START_REF] Guilbault | Enquête ECHO envoi -chargeurs -opérateurs de transport[END_REF], where a shipment is defined as a "quantity of freight that is made available at a given time in order to be transported during a single transportation operation from a given shipper to a given consignee", and consider it is as the relevant measurement unit with regard to the economics of freight transportation.

The main elements of the Economic Order Quantity model are now briefly recalled. Consider a firm sending a regular commodity flow of constant rate Q , from a given location to another by a given transportation mode. Freight transportation operations being discrete by nature, commodities are carried in shipments. We assume that all the shipments are of the same size s , and that each shipment is dispatched as soon as there are enough commodities at the origin. Then, the average origin stock level is 2 / s . It is the same at destination, provided that the commodity is sold or consumed at the rate it is produced. The average total inventory level is then s .

The shipper decides shipment size on the basis of all its costs. This concerns the freight rates paid to carriers, and the potential effort needed to carry out the transport operation, order costs, etc. These are assumed to consist of a fixed cost b independent of shipment size, and a variable cost s K. proportional to shipment size. Over a time period, s Q / shipments are sent by the shipper. Then, the total freight rates per time period amount to:

. s

Qb KQ + (1)
The willingness of the shipper to pay for a reduction in inventory is assumed proportional to the inventory level and to the time that elapses between the moment a commodity unit is produced and the moment it is sold, up to a coefficient a . The travel time is denoted by t . The total inventory is on average equal to Qt s + , so that the inventory cost per time period is ) ( Qt s a + . The total logistical cost function, denoted by g , is the sum of these two components:

. ) ( ) ( s Qb Q K at as s g + + + = (2)
The optimal shipment size is obtained by minimizing this convex cost function; it is the unique root of its first derivative:

( ) . * *) ( ' 2 s Qb a s g - = (3)
It depends only on a , Q , and b , and has the familiar "square-root" shape:

. * a Qb s =

* s does not depend on costs that are proportional to shipment size, such as the pipeline inventory cost ( atQ ) or the proportional component of the transportation cost ( K ). However, the total logistical cost does depend on such costs, and this has a direct effect and substantial effect on mode choice.

This model can easily be estimated by linear regression. Indeed, by taking the logarithm of both sides of Equation ( 4), we obtain:

a b Q s ln 2 1 ln 2 1 ln 2 1 * ln - + = ( 5
)
This equation is valid for each transport mode. As explained in the next section, Q and a are observed, and b ln ) 2 / 1 ( will be estimated as a mode-specific constant. From now on the optimal shipment size will simply be denoted by s , instead of * s .

DATA PRESENTATION AND MODEL SPECIFICATION

The ECHO database describes 10,462 shipments sent by some 3,000 shippers. These shipments are of all natures, sizes, origins and destinations (either the origin or destination is in France). All transportation modes are represented. Heavy modes, where the relative number of shipments is very low, were oversampled. There are much less observations than in other shipment databases, but each shipment is described with a large quantity of detail, which regards: the shipping and receiving firms; their relationships; the shipment itself; the transportation operation. In each of these categories, variables relating to economic, logistical or transport-related factors are available.

The equation that will be estimated is closely inspired by Equation [START_REF] Small | The economics of urban transportation[END_REF]. s ln is the dependent variable. Note that shipment size is measured in both weight and in volume in the ECHO database. The weight measurement is preferred in this paper; mainly because the volume measurement is much less accurate, insofar as it is given in m 3 , without decimals. Summary statistics of s and s ln can be found in Table 1. In Equation ( 5), the explanatory variables are Q , a and b . However, they are not observed directly in the ECHO database. Let us now examine how the variables available in the ECHO database can be used instead. For that purpose, variables relating to transportation demand and to transportation supply are addressed separately.

Transport demand related variables

Strictly speaking, the rate Q of the commodity flow between the shipper and the receiver is not available in the ECHO database. The closest available variable is the total commodity flow denoted by tot Q , without consideration of commodity type. This variable is available for 81.5% of the observations. Using tot Q instead of Q without care can lead to an underestimation of the influence of Q on s .

Nevertheless, there is a possibility that this bias may be limited. Indeed, in the ECHO survey, the shippers and receivers are identified at the level of facilities, or premises, i.e. physically well-defined components of firms. In addition, if a given facility is the origin (or destination) of a large number of distinct commodity flows, it is probable that these flows are being sent to (or coming from) distinct locations. In any case, it is probable that Q and tot Q often coincide, although no quantitative evidence can be drawn from the ECHO survey to support this statement. From now on tot Q will be used as a proxy for Q . The commodity value of time a is not directly available either. Fortunately, in this case, there is a good candidate to replace it. Indeed, for 64.5 % of shipments, the market value in euros (before tax) is specified. By combining this with the shipment weight, it is possible to calculate the value density of these shipments. The value density is denoted by dens a , and is measured in euro per kilogram.

Using dens a instead of a is a strong hypothesis. The commodity value of time, considered from the shipper's perspective, is certainly closely related to the value density. However, it is also expected to depend on many other parameters, such as the commodity's depreciation, the opportunity cost of capital for the shipper, the organization of the supply chain between the shipper and the receiver, and so on. Assume this can be represented by an interest rate applied to the value of the commodity. Then, due to the logarithmic specification of the model, this interest rate goes in the intercept of the model, and its variations in the residuals. Assuming this rate is not too strongly correlated with other exogenous variables in the model, this should not introduce a bias in the outcome of the estimation. Within the framework of this study, this assumption is made, and dens a is used as a proxy for a .

Transport supply related variables

The vast quantity of variables relating to transportation operations can be classified into two groups: technical variables, and freight rates. Technical variables describe the transportation technique and the availability and compatibility of other techniques. For example, live animals and chemicals cannot be carried in the same kind of trucks or rail cars. The technical variable used in the estimation is the main mode, that is to say the mode of the longest leg in the transportation operation (Table 2).

[ TABLE 2] Freight rates play a central role in the choice of shipment size. They are available in the ECHO database. However, the choice of shipment size actually depends on the generalized transportation cost from the perspective of shippers. The relation between freight rates and generalized transportation costs may be far from trivial. As a consequence, it has been decided not to use freight related variables in this study.. A simpler approach is preferred: on the assumption that b does not depend on origin-destination distance (given that b represents such costs as maneuvering, loading and unloading, and administrative operations that have no clear link with distance), it is sufficient to consider it as a modespecific constant (note that b can be assumed to depend on other variables, such as the number of legs in the transport operation and others; this is done in Section 6). Subsequently, the following specification is estimated:

1 + + + = mode mode mode dens a tot Q u X a Q s β β β ln ln . ( 6
)
where u is the error term and { } mode mode X are dummy variables indicating the transportation mode used.

ESTIMATION OF THE EOQ MODEL

The model is estimated using ordinary least square regression. The results are given in Table 3. The coefficients are highly significant. The R² coefficient is close to 0.8. Additional analyses of the residuals, available from the author or in the appendices of (45), do not invalidate the ordinary least square specification. Overall, specification [START_REF] Meyer | The economics of competition in transportation industries[END_REF] The absolute values of the mode-specific constants provide little information in themselves. However, their rankings are interesting. On average, the ranking of the modespecific constants should be the same as the shipment-independent component of generalized transportation costs b for each mode. Indeed, as shown in Table 3, the higher the mode capacity, the larger the mode-specific constant. It confirms the fact that the lower the unit transportation cost, the larger the fixed per-shipment cost, and the hierarchy which goes from waterway to rail, sea, combined transport, and then road and air transport, is intuitively correct. The apparently counter-intuitive gap between common carrier and private carrier is explained in the next section.

+ + + + + + + - = (7) Q β is close to 0.5,
[TABLE 3]

EXPLORATORY ESTIMATION OF AN EXTENDED EOQ MODEL

Other variables may be expected to influence shipment size. Some of them are available in the ECHO database, and their influence can be assessed by a small additional effort. In this section, the respective roles of the following variables are assessed: the origin-destination distance variable d , the number of agents physically or administratively intervening in the transportation operation interv N , the number of legs in the transportation operation trips N , and the organization of the transportation operation O (three possible organizations are distinguished: "isolated shipment", when the shipment is carried alone; "bundle", when the shipment is carried with other shipments; "part of a round", when the vehicle delivering or picking the shipment up reaches several distinct locations during the same movement). The latter two variables, when available, only concern the part of the transportation operation which has been carried out in the UE15 region. Summary statistics on these variables are given in Table 4. Variable O is summarized in Table 5.

[TABLE 4]

[TABLE 5]
In order to keep the estimation as simple as possible, and since there is no convincing microeconomic model relating b to these variables, a simple specification without any interactions is estimated:

u X X N N d X a Q s O round
+ + + + + + + - = 1 β β β β β β β β ln ln ln ln (8)
where O i X equals 1 if i O = , or otherwise 0; and u is the error term. Note that isolated β , the parameter corresponding to the case of isolated shipments, is set to zero by convention.

The model is once again estimated using ordinary least square regression. The results of the estimation are given in Table 6. A basic analysis of the residuals, available from the author or in [START_REF] Combes | The choice of shipment size in freight transport[END_REF], does not invalidate the ordinary least square hypotheses.

[ TABLE 6] All the coefficients are significant. Furthermore, as confirmed by the analysis of variance in Table 7, all the additional variables significantly improve the model. However, it remains true that the explanatory power of the model comes predominantly from tot Q ln and a ln , i.e. the core of the EOQ model. The estimated coefficients of tot Q ln and a ln remain reasonably close to (although not exactly the same as) the theoretical values. The hierarchy of the main transportation modes does not change substantially. It may be noted that the mode specific constants associated with private carriers and common carriers are not significantly different contrary to the previous model. The difference was due to the number of trips: the average number of trips when the main transportation mode is private carrier is 1.02, whereas it is 2.12 when the main transportation mode is common carrier.

Let us now briefly interpret the results obtained for the new variables. interv N has the expected positive effect: a larger number of agents intervening on the transportation chain seems to imply a larger b ln , thus a larger shipment size. What might be perceived as counter-intuitive is the negative sign of trips β . In fact, it confirms the economic rationale of hub-and-spokes transport networks: such networks are especially designed to handle small shipments efficiently. As a consequence, the transshipment cost is smaller in a hub-and-spokes transportation network than in more direct transportation organizations. Note that this is not in contradiction with the assessment that erv int β is positive; these two variables are scarcely related (their correlation is only 0.26).

[TABLE 7]

The influence of the transportation organization variable O on the shipment size is intuitive: shipments sent in bundles and in routes certainly share some fixed transshipment and handling costs, hence the negative signs in bundle β and round β . The positive influence of d ln on shipment size is perhaps the most difficult to understand. Normally, the distance should not influence the shipment size (subject to the use of a given mode). Here are three possible causes. The first is a frontier effect: shipments arriving from or going abroad entail more administrative work, which positively influences transportation cost per shipment. Second, unused vehicle capacity is more expensive to carriers on longer distances. It is possible that, through the freight rates, and particularly through surcharges and discounts, carriers induce shippers to send bigger shipments (46). Third, it may arise from the structure of logistics networks: if large flows travel from plants to regional distribution centers where they are then stocked for subsequent dispatch, there is a sort of logistical disconnect between the two operations (shipments are not targeted to a given retail center when they leave the plant, but only when they leave the regional distribution center), which results in the observation that bigger shipments move over longer distances.

CONCLUSION

In this paper, the Economic Order Quantity (EOQ) for shipment size choice in freight transportation is estimated using the French disaggregate shipment database ECHO. After a brief recap of the theoretical EOQ model, it is estimated using ordinary least squares. The estimation performs correctly, as about 80% of the variance in the sample is explained by the EOQ model. The theoretical model itself is clearly confirmed empirically.

From a microeconomic perspective, a number of conclusions can be drawn. First, inventory theory models, although (and perhaps because) they have been designed to optimize daily decision making for firms in well-defined logistical contexts, also prove valuable from the aggregate perspective of freight travel demand modeling. They offer theoretical insights into the behavior of shippers. In addition, it is reassuring to observe that a very simple model performs well at the aggregate scale of a large population of shippers carrying all kinds of commodities by all kinds of transportation modes. The intuition in [START_REF] Baumol | An inventory theoretic model of freight transport demand[END_REF] was relevant.

Second, shipment size is clearly dependent on transportation mode, and freight mode choice depends on shipment size as well. This paper brings an additional confirmation to the previous studies according to which shipment size choice and freight mode choice are simultaneous decisions.

Third, the significant role of a seemingly unimportant variable is highlighted in this paper: the rate of the commodity flow to which the shipment belongs. It plays an important role in determining shipment size and therefore most certainly also in determining mode choice. If this is true (further research is needed to confirm this statement, but it is consistent with the EOQ model, which fits very well in a mode choice framework) then -from the perspective of freight transportation demand -shippers should not only be characterized by their value of time and commodity type, but also by the rate of the commodity flows that they send to receivers.

Fourth, the availability of adequate data is crucial. It is widely agreed that logistics have a very powerful influence on freight transport, and need to be studied in detail. However, to do this, data are required, particularly on the logistical contexts in which transportation operations take place. Fortunately, the ECHO database contains such information, in addition to the characteristics of shipments themselves, and the related transportation operations, which are described with a lot of useful detail. It could be improved so that variables, such as the shipment volume, or the commodity flow between the shipper and the receiver concerning a particular commodity, be available.

This work offers many avenues for future research. The model used to explore shipment size choice could be improved in many ways; one of the most straightforward would be to account for the capacity constraints of vehicles, or to examine more closely how transportation prices, which are observed in the ECHO database, could be used in the estimation. Another natural direction for research would be to design a freight mode choice model, based on the EOQ model along the lines suggested in provided in [START_REF] Baumol | An inventory theoretic model of freight transport demand[END_REF] and [START_REF] De Jong | A micro-simulation model of shipment size and transport chain choice[END_REF]. This could lead to the development of a new type of freight transportation demand model, with an explicit (although still partial) representation of logistics. Other decisions, such as the choice of safety stock, could be introduced. This would entail significant efforts of both theoretical and empirical natures, as well as better surveys.

(46) Combes, F., and F. Leurent. On shipment size and freight tariffs: technical constraints and equilibrium price schedules, Proceedings of the 89 th Transportation Research Board Annual Meeting. Washington, USA, 2010.

LIST OF TABLES

  and a β relatively close to -0.5, as predicted by the theory. The basic EOQ model proves empirically effective in explaining the choice of shipment size by a large and heterogeneous population of shippers. It explains about 80% of the variance in the shipment size for which the explanatory variables are available (i.e. about 55% of the observations). The high importance of tot Q as an explanatory variable of shipment size is also confirmed. This tends to indicate that tot Q is an important explanatory determinant of mode choice in freight transport.

  seems adequate. The estimated model is:

	ln	s	0	.	50	ln	tot Q	0	.	44	ln	a	dens	1	.	05	X	commoncarr	ier
				1	.	46	X	rier privatecar			3	.	42	X	rail	2	.	09	X	combined
					4	.	37		X	waterway		2	.	89	X	sea	1	.	47	X	air

TABLE 1 Basic EOQ continuous variables summary statistics TABLE 2 Main transportation mode summary statistics TABLE 3 Estimation of the EOQ model TABLE 4 Extended EOQ continuous variables summary statistics TABLE 5 Shipment organization summary statistics TABLE 6 Estimation of the extended EOQ model TABLE 7 Analysis of variance in the extended EOQ model

 1

TABLE 1 Basic EOQ continuous variables summary statistics

 1

	Variable	Min.	Q1	Med.	Mean	Q3	Max	NAs
	1 (t)	0.00	0.05	0.65	19.58	7.8 10,800	0
	2343 (kt/y)	0.00	1.00	18.00	2.12 350.00 63,000	1,934
	56781 (k€/t)	0.00	1.07	4.56	59.37	20.00 10,400	3,715
	12313	-6.91	-3.00	-0.43	-0.65	2.05	9	0
	12323433	-6.90	0.18	2.89	3.02	5.86	13	1,934
	123567813	-0.94	6.97	8.43	8.49	9.90	16	3,715

TABLE 2 Main transportation mode summary statistics

 2

	9A	Number %
	Private carrier	1727 17
	Common carrier	6648 64
	Rail	224 2
	Combined transport	133 1
	Inland waterway	44 0
	Sea	825 8
	Air	859 8
	NAs	2 0

TABLE 3 Estimation of the EOQ model

 3

	Coefficients	Estimate	Std. Error	t-stat
	β	Q		0.50	0.01	73.27 ***
	β	a		-0.44	0.01 -37.57 ***
	β	common car	r rie	1.05	0.11	9.47 ***
	β	rrier private ca	1.46	0.11	12.87 ***
	β	rail		3.42	0.18	19.30 ***
	β	combined		2.09	0.20	10.31 ***
	β	waterway		4.37	0.33	13.05 ***
	β	sea		2.89	0.13	21.49 ***
	β	air		1.47	0.14	10.29 ***
	N			10,462
	NAs		4,741
	R 2		0.795
	Adjusted R 2	0.795

Significance levels: '.' at 10%; '*' at 5 %; '**' at 1%; '***' at 0.1 %

TABLE 4 Extended EOQ continuous variables summary statistics

 4

	Variable Min.	Q1	Med.	Mean	Q3	Max	NAs
	6 (km)	1.00	74.00 278.00 1,253.00 611.00 18,840.00	8
	42A63	0.00	4.30	5.63	5.44	6.42	9.84	8
	N	interv	0.00	1.00	2.00	2.76	3.00	12.00	720
	N	trips	1.00	2.00	2.00	2.06	3.00	8.00	720

TABLE 5 Shipment organization summary statistics

 5

	BA	Number %
	Isolated shipment	7647 73
	Part of a bundle	989 9
	Part of a round	1767 17
	NA's	59 1

TABLE 6 Estimation of the extended EOQ model

 6

	Coefficients	Estimate	Std. Error	t-stat
	β	Q		0.44	0.01	61.94 ***
	β	a		-0.43	0.01 -37.57 ***
	β	common car	r rie	0.95	0.12	7.57 ***
	β	rrier private ca	0.97	0.13	7.26 ***
	β	rail		2.78	0.19	14.56 ***
	β	combined		1.60	0.23	7.04 ***
	β	waterway		3.91	0.38	10.42 ***
	β	sea		1.59	0.19	8.57 ***
	β	air		0.34	0.18	1.85	.
				0.21	0.01	14.60 ***
	β	interv		0.16	0.02	8.92 ***
	β	trips		-0.40	0.02 -21.12 ***
	β	bundle		-0.62	0.07	-9.38 ***
	β	round		-0.69	0.06 -12.18 ***
	N			10,462
	NAs		5,134
	R 2		0.827
	Adjusted R 2	0.827

d β Significance levels: '.' at 10%; '*' at 5 %; '**' at 1%; '***' at 0.1 %

TABLE 7 Analysis of variance in the extended EOQ model

 7

	Variable Df		Sum Sq. Mean Sq. F-value
	ln	Q	tot	1 13,263.7 13,263.7	7452.8 ***
	ln	a	dens	1 28,915.7 28,915.7 16,247.5 ***
	M			7 1,563.8	223.4	125.5 ***
	ln	d			1	345.7	345.7	194.3 ***
	N	interv	1	95.8	95.8	53.8 ***
	N	trips	1	788.6	788.6	443.1 ***
	O				2	372.3	372.3	104.6 ***
	Residuals 5,314 9,457.3	1.8

Acknowledgements

The author would like to thank three anonymous referees for their helpful corrections and suggestions.


Jean-Baptiste Hoock
email: hoock@lri.fr

Jacques Bibai
email: jacques.bibai@gmail.com

Solving a Goal-planning task in the MASH project

Keywords: MASH, goal-planning, clustering, Monte-Carlo Tree Search

The MASH project is a collaborative platform with the aim to experiment different methods in an unknown environment of large size. The application is a goal-planning task in a 3D video game where runs are expensive. Moreover, there is no prior knowledge, the decisions have unknown semantics, observations on the environment are partial and of big size and accomplishing the task by taking random decisions always requires a very long run. So, solving this task is a big challenge. In this paper, we extend Monte-Carlo Tree Search, which has been proved very effective for applications in which simulating is easy and fast, to contexts in which there are only "real" expensive runs. This generic approach combines Clustering and Monte-Carlo Tree Search.

I. INTRODUCTION

Many specialized methods (e.g. in Data Mining [START_REF] Savaresi | Cluster selection in divisive clustering algorithms[END_REF] [START_REF] Defays | An efficient algorithm for a complete link method[END_REF]) have proven the efficiency in a lot of specific domains. Ensemble methods such as Adaboost [START_REF] Dubout | Tasting families of features for image classification[END_REF] try to take advantage of each specific method by building more collaborative approaches. In this context, the MASH project 1 was created. The aim of this project is to create new tools for the collaborative development of large families of feature extractors in order to start a new generation of learning software with great prior model complexity. Targeted applications of the project are (i) classical vision problems and (ii) goalplanning in a 3D video game and with a real robotic arm. For those applications, features are the result of one or several methods used to extract relevant information. The goal-planning application starts with a given initial state, a set of decisions and a set of methods. The state is the environment in which the avatar (for the 3D video-game) or the robot-arm are located. A decision is applied in order to go from a given state to another state. When methods are applied to a state, observations are generated. From these observations, the best decision is selected and then applied. A new state, a reward and an information about the final state are generated. A final state is the situation where the state doesn't change anymore whatever the applied decision. These processes are repeated from the initial state until one 1 http://mash-project.eu/wiki/index.php/About the MASH project. MASH means MAssive Sets of Heuristics final state is reached. The goal is to maximize the sum of all rewards given from an initial state to a final state. The problem is to find best decisions in order to reach the goal. Some algorithms such as Monte-Carlo Tree Search [START_REF] Kocsis | Bandit-based monte-carlo planning[END_REF], Nested Monte-Carlo Search [START_REF] Cazenave | Nested monte-carlo search[END_REF], Rollout Classification Policy Iteration [START_REF] Dulac-Arnold | Fast reinforcement learning with large action sets using errorcorrecting output codes for mdp factorization[END_REF] have been proposed. However, those algorithms can not be applied because either prior knowledge about decisions, environment and/or methods are crucial, either the complexity of the problem should be smaller (the number of features in the observation is not too large or the reward is informative...) either a model of the problem is required in order to simulate it. The MASH goal-planning problem is hard and there is no prior knowledge or model. This paper proposes a new generic algorithm in order to solve the goal-planning problem in the 3D video-game. The first part presents the MASH problem, the second part presents the algorithm and the third part presents experiments and discussions on 2 different problems and the last part is the conclusion.

II. THE MASH APPLICATION

We propose to solve a goal-planning problem in the 3D video-game. The environment is first described, then methods, observations, decisions and rewards. Finally, the best solution of this problem is given.

1) Environment: The environment is a square various size room with globally grey textures. The square room enclosed by 4 walls, contains an avatar and 2 flags. One of the two flags is red and the other flag is blue (See Fig. 1).

2) Methods: Methods are heuristics. More specifically, a heuristic is an image processing operator. The heuristic is applied on the view of the 3D avatar. 3 heuristics are used :

• red which indicates if a pixel (an element of the view of 3D avatar) is red or not • blue which indicates if a pixel is blue or not • identity which transforms the value of the pixel to a gray-scale value.

Each heuristic gives around 100,000 features. In the observation, the solver ignores what heuristic provides the feature and has no notion of heuristics.

3) Building an observation:

The cumulative number of features given by the 3 heuristics red, blue and identity are 3 × 100, 000 features but it's too huge. In order to reduce the size of observation space, 10,000 features are randomly selected over the 300,000 features. The observation given to the solver will be always made with these 10,000 features.

4) Decisions: There are 4 decisions :

• 0 (which means "go forward")

• 1 (which means "go backward")

• 2 (which means "turn right") • 3 (which means "turn left") Note that the solver doesn't know the meaning of decisions.

5) Rewards:

• Hit a wall: -1

• Touch the first time the blue flag: +5

• Touch (or hit) another time the blue flag: 0

• Hit the red flag without touched the blue flag: -5

• Hit the red flag after having touched the blue flag: +10

• Else: 0. This definition of rewards is not informative, because most often, the avatar moves without touching anything and receives a reward of 0.

6) Goal: A final state is reached if and only if once having touched the blue flag, the red flag is touched. The best cumulative reward is 15, by first touching the blue flag and then touching the red flag without hitting wall. We can distinguish 2 subgoals (i) touch the blue flag and then (ii) touch the red flag.

III. THE ALGORITHM

As it has been seen in the section II, there is no information that can help us to infer knowledge for solving more efficiently the problem. But we can simulate several scenarios in order to build a knowledge. After having defined Notations, 2 pieces of knowledge built by unsupervised methods are presented: (i) Categorization of Decisions and (ii) Clusterization of features. Then the Policy and the Algorithm, which uses built knowledge and learns how to reach the goal, are shown.

A. Notations

An observation (or a state) is denoted by s. From a given state s t , the state reached after the application of one decision on the state s t is denoted by s t+1 . The set of states is denoted by S, the set of features by F , the set of decisions by U , a feature by f and a decision by u. Let r t the reward and f inished t the Boolean variable which states if the state s t+1 is terminal (or final), M akeT ransition denotes the transition function (i.e. the application of the decision u t on the given state s t ) and returns the new state s t+1 , the boolean f inished t and the reward r t : (s t+1 , f inished t , r t ) = M akeT ransition(s t , u t ). When a transition is applied on a state s t , features of the

B. Learning a Categorization of the Decisions

At the beginning, the algorithm (or solver) has no idea of decisions. But by doing simulations, the impact on the state can be studied.

Formally, a decision u is :

• periodic if ∃k ∈ N * , ∀s t ∈ S s t+k = s t • stationary if ∀s t ∈ S ∃K ∈ N * , ∀j ∈ N * s t+K = s t+K+j • final if ∀s t ∈ S s t+1 is a final state.
with s t+x the state s t obtained after applying x times the decision u.

A decision u -is the inverse of a decision u + if s t = s t+2 with (s t+1 , , ) = M akeT ransition(s t , u + ) and (s t+2 , , ) = M akeT ransition(s t+1 , u -) ( denotes an ignored information.).

First, from a given state s t , we categorize a decision u by applying it several times. Let k the number of times the decision u has been applied.

1) A first case is to reach a final state by applying the first time (k = 1) the decision u. The decision is categorized as a final decision. 2) A second case is the decision has no more impact on the state (s t+k-1 ≈ s t+k ). The decision u is categorized as a stationary decision. 3) A third case is the initial state s t is seen a new time (s t+k ≈ s t ). The decision u is categorized as a periodic decision. 4) A fourth case is the initial state s t is seen a new time after applying the decision u one time (k = 1) and another decision u * = u one time, too. Decisions u and u * are categorized as inverse decisions. Pieces of information will be very useful for simulating more efficiently.

C. Learning a Clustering of the Features

The algorithm ignores what heuristic provides the feature. The idea is to clusterize relevant features that can help the algorithm to find goal. To solve this issue, we propose an algorithm of clustering (Fig. 2) called CluVo for finding correlated features. CluVo makes no assumption about the number of clusters and is completely unsupervised (Features are not labelled.). In the MASH problem, it corresponds to find, optimally, 2 clusters : "Features which come from the blue heuristic" and "Features which come from the red heuristic".

In the clustering, we distinguish 3 parts : (i) The generation of collection (or database) of lists of features by doing simulations -(ii) Clustering features -(iii) The Metric which is the correlation between 2 features.

1) Simulations for the generation of the collection: For the generation of the collection, following the number of iterations, simulations are done differently : For the first 3 iterations, decisions are chosen randomly; but from the 4 th iteration, decisions are chosen randomly or by vote (See Section III-D4). Given a cluster, the vote is very efficient for finding new features but needs several updates of (4) before to be used efficiently.

2) Clustering features: The approach of our algorithm is hierarchical. There exist 2 variants : (i) The divisive (or top-down) clustering [START_REF] Savaresi | Cluster selection in divisive clustering algorithms[END_REF] starts from one cluster and gradually splits clusters. (ii) The agglomerative (or bottomup) clustering [START_REF] Defays | An efficient algorithm for a complete link method[END_REF] starts from a set of small clusters and then aggregates them. Because of the system of vote (4), a bottom-up clustering is surely better. The shown algorithm is so agglomerative. Therefore, each feature corresponds to a cluster at the starting of the CluVo Algorithm.

3) Metric: Let ldb the collection of the list of features. Let f 1 and f 2 two features of the set of features obtained thanks to a uniform draw.

Let Db(f ) = {list ∈ ldb; f ∈ list} with f a feature and

Db(f 1, f 2) = Db(f 1) ∩ Db(f 2)
Let Card(X) = #X the number of elements of the set X.

The features f 1 and f 2 are significantly present (or correlated

) if #Db(f 1, f 2)) > 5 and 3 × #Db(f 1, f 2) >= max(#Db(f 1), #Db(f 2)).
The numbers 5 and 3 have been fixed after experiments.

D. Policy

The Policy π(s) is the function which gives the best decision to apply from the state s.

The best Policy should be able to go to the blue flag, first and then to go to the red flag. In order to realize this task, the Policy must be able to switch to the first subgoal (reach the blue flag) to the second subgoal (reach the red flag).

1) Memory for switching subgoals: A memory is necessary so as to memorize if the blue flag has been touched or not. The memory is modeled by a tree 2 of subgoals called ga. A node of the tree corresponds to a subgoal, an edge the decision of reaching a subgoal, the root means the top of the tree and a leaf is a node which has no next subgoal. A node contains (i) a rule of goal (ii) a list of features to activate (given by a cluster) with optionally, a number of features to activate (iii) the number of simulation, (iv) the cumulative reward of its last simulation and (iv) the average of cumulative rewards. 2 rules of subgoals are defined :

• GOAL AP P ROACH reach a state whose an amount quantity of features are activated • GOAL F IN ISHED the next decision should lead to a final state. To go from a subgoal gn 1 to another subgoal gn 2 , a first step (called IsReached(gn 1 , s t ) in Fig. 4) is to estimate if the state s t accomplishes the subgoal or not. If the answer is positive, a second step (called gn 2 = ChooseN extN ode(gn 1 ) in Fig. 4) is to choose from the current subgoal gn 1 the next subgoal to realize. The next chosen node is the child node of gn 1 which maximizes the number of simulations.

2) Useful Macro-Decisions: A Macro-Decision [START_REF] Mcgovern | Roles of macroactions in accelerating reinforcement learning[END_REF] is a decision repeated several times. A Macro-Decision is modeled by the couple (u, l).

• u is the decision • l is the number of times that the decision is repeated. A decision can be useless because the application of the decision on the state s t has no more impact (stationary decision), or the decision is the inverse of the last decision, or when the avatar has made a complete turn, it doesn't need to continue to turn (periodic decision). The set of useful decisions for a given state s t is denoted U ok .

3) Exploration for finding a subgoal: In the MASH application, the environment is partially observable. The environment must be explored. But because of the size of the environment, classical Random Search algorithm can not be applied. Macro-decisions is a tool that can help to explore more rapidly the environment [START_REF] Mcgovern | Roles of macroactions in accelerating reinforcement learning[END_REF]. To touch a flag, a first step consists to find it.

To find the flag, some decisions are more useful than others (e.g. the decision to turn is good). During the exploration, when there is no activated features in the state s t and then a given decision activates features in the state s t+1 , this decision can be considered as a good one. These decisions can be found statistically after a lot of explorations.

In order to choose good decisions for finding a subgoal, 2 variables are defined. Let F clus a set of features

disc(F clus, u) = #{s t , F a (t) == ∅∧F a (t+1)∩F clus = ∅} (1) disc(u) = #{s t , F a (t) == ∅ ∧ F a (t + 1) = ∅}. (2)
The formulas can be updated at each simulation after a M akeT ransition.

Then for a given node gn, the decision for exploration is given by • Let obj the set of features given by the node gn.

• one time over 1 + u∈U disc(u), return a random decision.

• with a probability of disc(obj, u)/disc(u), if u ∈ U ok then return the decision u • else return a random decision When features have been activated, decisions can be now chosen by using these activated features.

4) Vote for reaching a subgoal: Given a set of features, a set of decisions and from an initial state with some active features, a sequence of decisions which leads to the goal is searched. A decision can activate new features or deactivate features. Here, the goal is to activate a maximal number of features of the set. In order to do this, the idea is that each of active features votes for a decision which increases the number of active features of the set.

For each couple (f, u) ∈ F × U is associated a score score(f, u). A vote is defined by

V ote(f ) = Argmax u∈U score(f, u).
(

) Let F v (u) = {f ∈ F a (t) ∩ F Clus, u == V ote(f ) 3
} the set of active features of the state s t belonging to a cluster (or set of correlated features) F Clus and which vote for the decision u.

Then, the decision by vote is given by Argmax u∈U ok #(F v (u)).

We can see like this : For each decision u, we count the number of active features which vote for the decision u. The decision which has received the most of votes is selected.

The function score gives a cumulative number of features of one cluster which have been activated from a sample of states. More precisely, Let (s t , u, s t+1 ) the state s, the decision u and the state s t+1 obtained immediately after applying the decision u on the state s.

Let F a|F clus (t) the set of active features of the state s t belonging to the cluster F clus (F a|F clus = F a (t) ∩ F clus). Let f a ∈ F clus an active feature of the state s. Then score(f a , u) = st (#F a|F clus (t + 1) -#F a|F clus (t)).

(4)

The signification of the formula is to activate more features in the state s t+1 than activated features of the state s t .

Note that f a ∈ F a (t) but maybe, f a / ∈ F a (t + 1). In this case, the feature f a has been deactivated when the decision u has been applied on the state s t .

The function score can be updated at any moment when a decision is applied.

E. Learning a Sequence of Goals

The categorization of decisions and the Clustering of features are unsupervised methods, the information of the reward has not been used. Now, for finding the goal, the reward is used and we build the tree ga. The proposed algorithm (Fig. 5) is a Monte-Carlo Tree Search [START_REF] Kocsis | Bandit-based monte-carlo planning[END_REF] [2] on subgoals. The algorithm is called GMCTS (Goal Monte-Carlo Tree Search).

Function Gmcts(maxSimulation) Let ga a static variable containing a tree of subgoals Let root the root of the tree ga. for sim in [[1; maxSimulation]] do Let s initialized to the initial state.

(rcum, f inished, gn) = (0, f alse, root) while f inished = true do gn = ChooseN extGoal(gn) reached = f alse while reached = true ∧ f inished = true do u = GetDecision(s, gn) l = GetM acroDecisionLength(s, u, gn) for j in 1..l do (s, f inished, r) = M akeT ransition(s, u) rcum = r + rcum reached = IsReached(gn, s) if reached ∨ f inished then break end if end for{this loop is the application of the Macro-decision (u, l) on the state s} end while end while U pdateT ree(rcum, gn) {Pieces of information in a node which are to update are : (i) the number of simulations (ii) the average of the cumulative rewards and (iii) the cumulative reward. All nodes visited during the simulation (going from the node gn to root) are updated.} end for Figure 5. Monte-Carlo Tree Search for searching goal. The argument maxSimulation is a parameter of the GMCTS algorithm and denotes the number of simulations which will be done to build the tree ga. Macrodecisions are used in the GMCTS algorithm. The policy of the GMCTS is given by the Fig. 6 Given the node gn, the function ChooseN ewGoal creates a new child node of gn and returns it, whereas the function ChooseBestGoal returns the child node of gn which maximizes either the average of cumulative rewards (mode 1) or the cumulative reward of its last simulation (mode 2).

IV. EXPERIMENTS

A. The optimization of a draw of letters

For the genericity, the Algorithm has been applied on another application : optimization of a draw of letters. The longest word is a game where a draw of 10 letters is made and the goal is to produce the longest correct word according to a dictionary of around 300,000 French words 3 . But, in order to have interesting draws, letters are not drawn equiprobably. The goal is to find the best distribution of letters in order to have draws with long correct words. 3 http://www.pallier.org/ressources/dicofr/dicofr.html Results of GMCTS on LD(20,100,10). Run 1 and Run 2 correspond to 2 runs. The mean cumulative reward corresponds to the average length of longest words on draws of 10 letters. State: The state is a vector of size N bL × 26; each of the 26 groups contains N bL features. Each feature is 0 or 1.

The state describes a distribution on words as follows: we randomly draw 10 letters, and we start by randomly draw a consonant -each consonant has probability proportional to the number of 1 in its features. Then, for each letter, we switch from vowel to consonant or from consonant to vowel with probability 95%; and letters are drawn among consonants or vowels with probability proportional to the number of features at 1. The initial distribution contains each letter one times.

Decisions: 27 decisions are possible. Each decision (except the decision "not add") consists in adding one letter; the feature which is activated among the N bL corresponding features is randomly drawn.

Reward: The reward is 0 if s is not terminal. When s is terminal, the reward is computed as the average length of longest words through N bDraws draws.

Terminal State: Let M al the maximal number of letters we can add in all. A state is terminal if either M al letters have been added either the decision of "not add" has been taken.

For different variants of this problem, the letter's draw application will be denoted LD(N bL, M al, N bDraws). A difficulty of this testbed is that the reward is stochastic. Like in the Mash application, we work in a setting with no prior knowledge.

B. Results

For both applications, all tests have been made with 100 evaluations and all experiments have been exactly made with the same algorithm. No parameter has been modified.

1) Results on draw of letters: Results of clustering is given in Table I. Results of GMCTS are shown in Figures 7 and8 in which the number of iterations is the number of simulations of the GMCTS Algorithm. The mean cumulative reward corresponds to the average length of longest words.

2) Results on Mash application: Results of the Clustering are : Over 12 runs, the number of runs with a bad correlation of features has been of 9. For the 3 other runs, 2 runs needed 5 iterations, one needed 8 iterations. At the end of the Clustering, 2 clusters have been well found.

Results of GMCTS are presented in the Table II.

3) Discussions: Letter's draw optimization: The clustering works well. Although in LD(100, 2000, 10) application, the size of the state space is 5 times larger and the maximal length of a simulation is 10 times bigger than in LD(20, 100, 10), the clustering works always well and the GMCTS Algorithm doesn't need more simulations to converge in average. In a similar application, the french TV game "Des Chiffres et Des Lettres", the longest word has in average a length of 8.12 letters (Average calculated over 275 draws). In Figures 7 and8, best runs have very good results; their score (given by the mean cumulative reward) reaches or exceeds a little 8.12. The score in average (around 7.5) is good, too.

Mash Application: The results are slow but stable. With 12 restarts on the clustering stage (2 hours each), we get 3 successes (which were detected as successes by heuristics without further testing); then, GMCTS could be successful on each of these 3 runs (8 hours each run). So on a parallel machine all this could be done in 10 hours, or 48 hours on a sequential machine. The success rate of CluVo is 25% (tested on 12 runs, and success can be detected on the fly) and the success rate of GMCTS is 100% (i.e. the goal "touch the blue flag then touch the red flag" is successfully found).

V. CONCLUSION

In this paper, we propose a generic algorithm CluVo+GMCTS which solves a very complex task in an unknown environment of large size. In this environment, the semantics of decisions are unknown, there is no prior knowledge and pure random search gives very rarely a good solution. Our approach combines Clustering and Monte-Carlo Tree Search. The Clustering method tries to learn some knowledge and from those results, GMCTS tries to find a strategy for reaching the goal. The genericity of the Algorithm has been shown on 2 applications. In further work, we will compare GMCTS with other algorithms such as Nested Monte-Carlo Search adapted to goal space.

VI. ACKNOWLEDGMENT:

Figure 1 .

 1 Figure 1. Top: 3 different views of the avatar. Bottom: Images resulting of the application of the blue heuristic. Left: The blue flag is not seen, no features are activated. Center: The blue flag is far, few features are activated. Right: The blue flag is closer, more features are activated state s t are activated or deactivated. Given a state s t , F a (t) denotes the set of active features of this state.

Figure 2 .

 2 Figure 2. the complete Clustering (generation of the collection + clustering).

Figure 3 .Figure 4 .

 34 Figure 3. The Clustering phase. SelectT woF eatures(lclus, ldb) selects 2 features by using either the list of clusters lclus or the collection of the list of features ldb. M erge(clus1, clus2, lclus) creates a new cluster into the set lclus by merging clus1 and clus2 and then removes them from the set lclus.Function π(s) Let ga a static variable containing a tree of goals Let gn the current node and root the root node of ga Let u a static variable containing the decision Let l a static variable containing the remaining number of times that the decision u will be repeated if gn is root ∨ IsReached(gn, s) then gn = ChooseN extN ode(gn) l = 0 end if if l == 0 then Build a new Macro-Decision (u, l) • u = GetDecision(s, gn)

Figure 6 .

 6 Figure 6. The Policy of GMCTS.

Figure 7 .

 7 Figure 7.Results of GMCTS on LD(20,100,10). Run 1 and Run 2 correspond to 2 runs. The mean cumulative reward corresponds to the average length of longest words on draws of 10 letters.

Figure 8 .

 8 Figure 8. Results of GMCTS on LD(100,1000,10). The errorline of the average curve is the standard deviation computed on all remaining runs. The errorline of the curve "One Run" is the standard deviation computed through 100 evaluations.

  Function CluV o() Let clusters a static variable containing all clusters of features. Initialization: each cluster is composed of one feature. for iteration in 1.. + ∞ do ldb = GenerateN ewDatabase(iteration) nbClus = size(clusters) clusters = Clusterize(ldb, clusters)

	if nbClus = size(clusters) then
	nbClus = size(clusters)
	else
	break
	end if
	end for

  Function Clusterize(ldb, lclus) Let P C1 a static parameter fixed to 10 4 Let P C2 a static parameter fixed to 10 7 Let nbClus the number of clusters initialized to the size of ldb nbT ry = min(max(P C1, nbClus × nbClus), P C2) while nbT ry > 0 ∧ nbClus > 1 do Decrement nbT ry (f 1, f 2) = SelectT woF eatures(lclus, ldb) Let clus1 the cluster of features to which the feature f 1 belongs Let clus2 the cluster of features to which the feature f 2 belongs if clus1 = clus2 ∧ f 1 and f 2 are correlated (See Section III-C3) then

	M erge(clus1, clus2, lclus)
	nbClus = size(lclus)
	nbT ry = min(max(P C1, nbClus × nbClus), P C2)
	end if
	end while
	Return lclus

Table I RESULTS

 I OF CLUSTERING ON LETTER'S DRAW APPLICATION. THE ROW ITERATIONS IS THE NUMBER OF ITERATIONS OF CLUSTERING. THE ROW ERROR IS THE NUMBER OF RUNS WITH BAD CORRELATION(S). THE ROW INCOMPLETE IS THE NUMBER OF RUNS WHERE THE CORRECT NUMBER OF CLUSTERS HAS NOT BEEN FOUND. THE ROW SUCCESS IS THE NUMBER OF RUNS WITH COMPLETE CLUSTERING WITHOUT BAD CORRELATION.

	Application	LD(20,200,10)	LD(100,2000,10)
	Nb Runs	43	35
	Iterations	between 3 and 4 3 (except one run with 4)
	Error	2	0
	Incomplete	1	0
	Success	40/43	35/35
		Table II

RESULTS

  OF GMCTS ON MASH APPLICATION. THE COLUMN ITERATION IS THE NUMBER OF SIMULATION OF GMCTS. AVGCR MEANS THE AVERAGE OF CUMULATIVE REWARD. THE OPTIMAL CUMULATIVE REWARD IS 15. THE LAST COLUMN SUCCESS IS THE PERCENTAGE OF EVALUATIONS WHERE THE GOAL "TOUCH THE BLUE FLAG AND THEN THE RED FLAG" HAS BEEN ACCOMPLISHED (THROUGH 100

	Run	Iteration	AvgCR	Success
	1	1	-13.65
		10	4.82
		100	9.67
	2	1	8.87
		10	10.62
		100	10.2
	3	1	12.09 ± 1.06	81%
		10	12.64 ± 1.03	86%
		100	11.38 ± 1.21	77%

a definition can be found in http://en.wikipedia.org/wiki/Tree (graph theory)

The authors would like to thank the whole Mash team who has permitted us to make these experiments and who has contributed to develop the algorithm through many discussions and Olivier and Fabien Teytaud for the reading of this article. Experiments on the application of the draw of letters presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr).


Cédric Simenel
email: cedric.simenel@anu.edu.au

Aditya Wakhle

Benoît Avez

Quantum microscopic approach to low-energy heavy ion collisions

The Time-dependent Hartree-Fock (TDHF) theory is applied to the study of heavy ion collisions at energies around the Coulomb barrier. The competition between fusion and nucleon transfer mechanisms is investigated. For intermediate mass systems such as 16 O+ 208 Pb, proton transfer favors fusion by reducing the Coulomb repulsion. A comparison with sub-barrier transfer experimental data shows that pairing correlations are playing an important role in enhancing proton pair transfer. For heavier and more symmetric systems, a fusion hindrance is observed due to the dominance of the quasi-fission process. Typical quasi-fission time of few zeptoseconds are obtained. Actinide collisions are also investigated both within the TDHF approach and with the Ballian-Vénéroni prescription for fluctuation and correlation of one-body observables. The possible formation of new heavy neutron-rich nuclei in actinide collisions is discussed.

Introduction

Microscopic description of nuclear reactions at low energy where quantum effects play a significant role is an important challenge of nuclear physics. The interplay between nuclear structure and reaction mechanisms is crucial at energies around the Coulomb barrier generated by the competition between Coulomb and nuclear interactions. For instance, this has been well established for fusion reactions (for a review, see, e.g., Ref. [1]).

A good understanding of these effects is crucial for the quest of superheavy elements (SHE) [2]. Indeed, different structures of the collision partners may change SHE production cross sections by several orders of magnitudes. For instance, collisions of quasi-symmetric heavy systems with typical charge products Z 1 Z 2 ≥ 1600 exhibit a fusion hindrance when compared to systematics and models established for lighter systems [3]. The dynamical shape evolution of the dinuclear system and the transfer between the fragments are expected to play a significant role on this fusion hindrance.

In the case of the collision of actinide nuclei, the Coulomb repulsion is so strong that fusion cannot occur. However, a dinuclear system may survive up to few zeptoseconds, forming the heaviest nuclear systems available on Earth. Multi-nucleon transfer between actinides may provide an alternate way to produce more neutron rich heavy and super-heavy elements. The complex dynamics of such systems with up to ∼ 500 nucleons in interaction is very challenging and has been investigated with different theoretical approaches recently [4,5,6,[START_REF] Zagrebaev | Journal of Physics: Conference Series (JPCS)[END_REF][START_REF] Maruyama | [END_REF]9,10,11,12].

The time-dependent Hartree-Fock (TDHF) theory [13] presents the advantage of treating both structure and reaction dynamics within the same formalism. The TDHF theory provides a self-consistent mean-field description of the many-body dynamics. Early TDHF calculations in nuclear physics used various symmetries and simplified Skyrme interactions [14] to reduce computational times [15,16,17]. Recent increase of computational power allowed realistic TDHF calculations in 3 dimensions [18,19,20,21,22] with full Skyrme energy density functional (EDF) including spin-orbit terms [14,23]. Several examples of recent applications of TDHF to nuclear dynamics were presented at the NN2012 conference [START_REF] Itagaki | Journal of Physics: Conference Series[END_REF][START_REF] Iwata | Journal of Physics: Conference Series (JPCS)[END_REF][START_REF] Oberacker | Journal of Physics: Conference Series (JPCS)[END_REF][START_REF] De La Mota | Journal of Physics: Conference Series[END_REF].

In Section 2, we recall the TDHF formalism and provide some details on the calculations. In Section 3, we study the interplay between fusion and transfer in the 16 O+ 208 Pb system. Section 3 is dedicated to the description of sub-barrier transfer in this reaction. In Section 5, we investigate the fusion hindrance of heavy quasi-symmetric systems and the quasi-fission process. Finally, we study the collision dynamics of actinide nuclei in Section 6 before to conclude.

The Time-Dependent Hartree-Fock approach

In nuclear physics, the TDHF theory is applied with a Skyrme EDF modeling nuclear interactions between nucleons [14]. The EDF is the only phenomenological ingredient of the model, as it has been adjusted on nuclear structure properties like infinite nuclear matter and radii and masses of few doubly magic nuclei [23]. The main approximation of the theory is to constrain the manybody wave function to be an antisymmetrized independent particle state at any time. The latter ensures an exact treatment of the Pauli principle during time evolution. Though TDHF does not include two-body collision term, it is expected to treat correctly one-body dissipation which is known to drive low energy reaction mechanisms as Pauli blocking prevents nucleon-nucleon collisions.

The main advantage of TDHF is that it treats static properties and dynamics of nuclei within the same formalism and the same EDF. The initial state is obtained through static HF calculations which give a good approximation of nuclear binding energies and deformations [START_REF] Vautherin | [END_REF]29]. Another important advantage of TDHF for near-barrier reaction studies is that it contains all types of couplings between the relative motion and internal degrees of freedom. However, TDHF gives only classical trajectories for the time-evolution and expectation values of one-body observables. In particular, it does not include tunneling of the many-body wave function and may underestimate fluctuations of one-body observables [30,31]. The latter fluctuations can be estimated in the time-dependent RPA limit using a prescription proposed by Balian and Vénéroni [32].

Inclusion of pairing correlations responsible for superfluidity in nuclei have been done recently to study vibrations in nuclei [33,34,35,36]. However, realistic applications to heavy ion collisions are not yet achieved and are beyond the scope of this work.

Formalism

The TDHF equation can be written as a Liouville-Von Neumann equation

ih ∂ ∂t ρ = [h[ρ], ρ] (1)
where ρ is the one body density matrix associated to the total independent particle state with elements

ρ(rsq, r ′ s ′ q ′ ) = A 1 +A 2 i=1 ϕ i (rsq) ϕ * i (r ′ s ′ q ′ ). (2)
The sum runs over all occupied single particle wave functions ϕ i and r, s and q denote the position, spin and isospin of the nucleon, respectively. The Hartree-Fock single particle Hamiltonian h[ρ] is related to the EDF, noted E[ρ], by its first derivative

h[ρ](rsq, r ′ s ′ q ′ ) = δE[ρ] δρ(r ′ s ′ q ′ , rsq) . (3)

Practical aspects

A TDHF calculation of two colliding nuclei is performed assuming that the two collision partners are initially at a distance D 0 in their HF ground state. This distance has to be large enough to account for Coulomb excitation in the entrance channel (polarization, vibration, rotation...). Typical initial distances are within the range D 0 ∼ 30 -50 fm. We assume that before the initial time, the nuclei followed a Rutherford trajectory. This assumption determines the initial velocities in the center of mass frame. The collision partners are put into motion using a Galilean boost [37] which is applied on each nucleus at the first iteration. We use the tdhf3d code built by P. Bonche and coworkers with the SLy4d Skyrme parametrization [18]. This code has a plane of symmetry (the collision plane). It uses the Skyrme energy functional expressed in Eq. (A.2) of Ref. [38] where the tensor coupling between spin and gradient has been neglected. The lattice spacing is ∆x = 0.8 fm and the time step is ∆t = 1.5 × 10 -24 s. More details on modern TDHF calculations to nuclear dynamics can be found in Refs. [START_REF] Simenel | Quantum Many-Body Dynamics: Applications to Nuclear Reactions[END_REF][START_REF] Simenel | [END_REF].

Fusion with medium mass systems

A natural application of the TDHF approach in nuclear physics is to study fusion reactions at and above the barrier. Fusion occurs by transferring relative motion into internal excitation via one-body mechanisms well treated by the TDHF approach. As a result, early TDHF codes have been successfully applied to describe above-barrier fusion reaction in light systems [41]. In this section, we investigate the fusion process in the 16 O+ 208 Pb system.

A reference nucleus-nucleus potential could be obtained in the frozen approximation with HF (or HFB) densities [42], where the collision partners are assumed to keep their ground-state density during the approach. The frozen potential can be computed with the same Skyrme EDF as in the TDHF calculations by translating the nuclei in their HF state. A comparison between TDHF and frozen fusion barriers allows to identify the role of dynamical effects, which are included in TDHF but absent in the frozen approach.

Writing the HF energy E[ρ] as an integral of an energy density H[ρ(r)], i.e.,

E[ρ] = dr H[ρ(r)], (4)
we get the expression for the frozen potential

V (R) = dr H[ρ 1 (r) + ρ 2 (r -R)] -E[ρ 1 ] -E[ρ 2 ], ( 5
)
where R is the distance between the centers of mass of the nuclei, and ρ 1,2 are the densities of their HF ground-state. Eq. ( 5) neglects the Pauli principle between the nucleons of one nucleus and the nucleons of the other one. As a result, it is valid only for small overlaps between the two collision partners. However, for light and intermediate mass systems, the barrier radius is large enough to justify this approximation. The HF-frozen potential for the 16 O+ 208 Pb system is shown in Fig. 1. The potential obtained from the Wong formula [43] is also presented for comparison. The height of the barrier obtained We now investigate the possible role of dynamical effects on the fusion barrier with TDHF calculations. The TDHF fusion barrier is identified as the capture threshold for central collisions, above which a compound system is formed and below which two fragments are emitted. Due to the finite time of the TDHF evolutions, one has to define a maximum computational time (∼ 10 3 fm/c in the present case) above which the final configuration (i.e., one compound system or two fragments) is assumed to be reached.

Figure 3 shows the evolution of the relative distance between fragment centers of mass in central 16 O+ 208 Pb collisions at different energies around the barrier. We identify two sets of trajectories associated to fusion for E ≥ 74. [START_REF] Guo | proceeding of the NSRT12 conference[END_REF] MeV and to re-separation of the fragments for E ≤ 74.44 MeV. The resulting fusion barrier predicted by these TDHF calculations is, then, V T DHF B = 74.445 ± 0.005 MeV. As a result, the dynamical effects included in TDHF calculations lower the barrier by ∼ 1.5 MeV for this system as compared to the HF-frozen approximation. Other methods based on a macroscopic reduction of the mean-field dynamics, namely the dissipative-dynamics TDHF [44] and the density-constrained TDHF [48], also find similar results with an energy dependence to the barrier heights ranging from 74.5 at low energies to 76 MeV at higher energies where the frozen approach is more valid.

The density evolutions close to the fusion threshold, i.e., at E c.m. = 74.44 and 74.45 MeV, are plotted in Fig. 4. In both cases, a "di-nuclear" system is formed during a relatively long time (∼ 500 fm/c) before either re-separation or fusion. Nucleon transfer is expected to occur in this di-nuclear system. An evidence of this transfer is given by integrating the proton and nucleon densities at the final time when the two fragments re-separate at E c.m. = 74.44 MeV. As a result, the outgoing channel is, in average, 14 C+ 210 Po, showing the importance of proton transfer from the light to the heavy partner, in good agreement with experiment [49,50]. This proton transfer effectively lowers the barrier by decreasing Z 1 Z 2 and, then, the Coulomb repulsion. As a result, the coupling to transfer channels increase the fusion probability for this system. Transfer reactions in the 16 O+ 208 Pb system are discussed in more details in the next section. Note that low-lying collective vibrations are also known to affect the fusion barrier distribution [47].

Transfer reactions

Transfer reactions in heavy-ion collisions have been investigated within the TDHF framework recently [51,52,53,[START_REF] Simenel | Quantum Many-Body Dynamics: Applications to Nuclear Reactions[END_REF]54,55]. In fact, the previous study showed the importance of the interplay between fusion and transfer reactions in the 16 O+ 208 Pb system around the barrier. For instance, just below the barrier, the light fragment in the exit channel (see top of Fig. 4) is, in average, a 14 C nucleus, indicating that the transfer of two protons is the dominant channel at this energy. We now discuss this sub-barrier transfer reaction in more details.

In the previous section, a signature for transfer mechanism was obtained from the change of average numbers of nucleons in the fragments. Using the quantum nature of the TDHF theory, it is possible to compute the probability for a specific outgoing channel, i.e., the probability to find a given number of neutrons and protons in the fragments [56,54]. This is possible, e.g., thanks to a particle number projection technique applied on the outgoing fragments [54].

Let us consider a final time when the two fragments are located on both side of the x = 0 plane and define a particle number projector onto N protons or neutrons in the x > 0 region

PR (N ) = 1 2π 2π 0 dθ e iθ( NR -N ) , (6)
where NR = s dr â † (rs) â(rs) Θ(x) [START_REF] Zagrebaev | Journal of Physics: Conference Series (JPCS)[END_REF] counts the number of particles in the x > 0 region (Θ(x) = 1 if x > 0 and 0 elsewhere). Isospin is omitted to simplify the notation. The projector defined in Eq. ( 6) can be used to compute

where |φ R (θ) = e iθ NR |φ . Note that |φ R (θ) is an independent particle state and, then, the last term in Eq. ( 8) is the determinant of the matrix of the occupied single particle state overlaps [54]:

φ|φ R (θ) = det(F ) (9)
with

F ij = s dr ϕ s i * (r)ϕ s j (r)e iθΘ(x) . (10)
The integral in Eq. ( 8) is discretized using θ n = 2πn/M with the integer n = 1 • • • M . Choosing M = 300 ensures numerical convergence for the 16 O+ 208 Pb system. Fig. 5 shows the resulting transfer probabilities at E c.m. = 74.44 MeV (left) and at E c.m. = 65 MeV (right). In agreement with the results presented in section 3, the most probable channel is Z = 6 and N = 8 at the barrier. However, lowering the energy reduces the transfer probabilities and the main channel is Z = N = 8 well below the barrier, corresponding to inelastic and elastic channels.

To compare with experimental data on the 16 O+ 208 Pb reaction, we plot in Fig. 6 the transfer probabilities from TDHF as a function of the distance of closest approach R min between the collision partners [57], assuming a Rutherford trajectory [START_REF] Broglia | Heavy Ion Reactions: Lecture Notes the Elementary Processes (Frontiers in Physics)[END_REF]:

R min = Z 1 Z 2 e 2 [1 + cosec(θ c.m. /2)]/2E c.m. ( 11
)
where θ c.m. is the center of mass scattering angle. Recent data from Ref. [50] are shown in Fig. 6 for sub-barrier one-and two-proton transfer channels. We see that TDHF overestimates the one-proton transfer probabilities and underestimates the strength of the two-proton transfer channel. This discrepancy is interpreted as an effect of pairing interactions [54,50]. Indeed, paired nucleons are expected to be transferred as a pair, increasing (resp. decreasing) the two-nucleon (single-nucleon) transfer probability. For R min > 13 fm, however, the TDHF calculations reproduces reasonably well the sum of one and two-proton transfer channels. For R min < 13 fm, sub-barrier fusion, not included in the TDHF calculations, reduces transfer probabilities [54,50]. These studies emphasize the role of pairing interactions in heavy-ion collisions. The recent inclusion of pairing interactions in 3-dimensional microscopic codes [34,35,36] gives hope in our ability to describe such data more realistically in a near future.

Fusion and quasi-fission in heavy systems

To a reasonably good approximation, the fusion barrier for light and intermediate mass systems is determined by the frozen barrier. This approximation fails, however, for heavy systems with typical charge products Z 1 Z 2 ≥ 1600 which are known to exhibit fusion hindrance [3]. Above this threshold, an extra-push energy is usually needed for the system to fuse [START_REF] Swiatecki | [END_REF]. In fact, at the energy of the frozen barrier, heavy systems are more likely to encounter quasi-fission, i.e., a re-separation in two fragments after a possible mass exchange.

We now investigate the reaction mechanism in heavy systems with possible fusion hindrance. We first illustrate the fusion hindrance with TDHF calculations of fusion thresholds. Then, we investigate the quasi-fission process.

TDHF calculations of fusion hindrance

Let us first consider the 90 Zr+ 124 Sn system which has a charge product Z 1 Z 2 = 2000, and, then, is expected to exhibit a fusion hindrance. Indeed, the proximity model [60] predicts a barrier for this system V prox. ≃ 215 MeV, while TDHF calculations predicts that the system encounters a fast re-separation at this energy, as shown in Fig. 7 [61]. In the same figure, we observe a long contact time at E c.m. = 240 MeV, which is interpreted as a capture trajectory leading to fusion. The additional energy needed to fuse is then ∼ 25 MeV, which is higher than the extra-push model [START_REF] Swiatecki | [END_REF] prediction E X-push ≃ 14.8 MeV.

TDHF calculations of fusion hindrance have been performed for several systems, as shown in Fig. 8, where the TDHF fusion thresholds are compared with the interaction barriers predicted by the proximity model [60] and with the extra-push model [START_REF] Swiatecki | [END_REF]. We observe an increase of . TDHF fusion thresholds for several heavy systems are compared with the proximity barrier [60] and with results from the extra-push model [START_REF] Swiatecki | [END_REF].

the fusion threshold with TDHF as compared to the proximity model which assumes frozen reactants. This indicates that dynamical effects are playing an important role in the reaction by hindering fusion. The order of magnitude of the additional energy needed to fuse is similar to the one predicted with the extra-push model.

The quasi-fission process

The quasi-fission mechanism is a fast re-separation of the fragments, with usually a partial equilibration of their mass. Typical quasi-fission times are of the order of few zeptoseconds [62,63,[START_REF] Simenel | [END_REF]. Quasi-fission becomes dominant in heavy systems and is mostly responsible for the fusion hindrance discussed in the previous section. As an example, Fig. 9 shows the density profiles for the 90 Zr+ 124 Sn head-on collision at E c.m. = 235 MeV. We observe that the two fragments are in contact during ∼ 5 zs. A dinuclear system is then formed before re-separation in two fission-like fragments.

Extensive calculations are ongoing on the 40 Ca+ 238 U system to compare with recent measurements performed at the Australian National University [64]. Examples of density evolutions are shown in Fig. 10 for two different initial conditions. In both cases, a quasifission process is obtained. In particular, we observe an important multi-nucleon transfer from the heavy fragment toward the light one. Fig. 11 presents final fragment masses (top) and quasi-fission times (bottom) for two different orientations of the 238 U. Comparing these two figures, we observe that the mass equilibration (i.e., the formation of two fragments with symmetric masses) is not complete and varies with the life-time of the dinuclear system, i.e., the longer the contact, the larger the mass transfer. We also see that all the calculations with the 238 U deformation axis aligned with the collision axis lead to a quasi-fission with partial mass equilibration and quasi-fission times smaller than 10 zs. Shell effects may affect the final outcome of the reaction by favouring the production of fragments in the 208 Pb region. In particular, this orientation never leads to fusion, while the other orientation produces long contact times above the barrier which may be associated to fusion. Calculations of non-central 40 Ca+ 238 U collisions are ongoing in order to compare with experimental data.

Actinide collisions

Collisions of actinides form, during few zs, the heaviest nuclear systems available on Earth. In one hand, such systems are interesting to study the stability of the QED vacuum under strong electric fields [65,66,11]. In the other hand, they might be used to form neutron-rich heavy and super-heavy elements via multi-nucleon transfer reactions [6,[START_REF] Zagrebaev | Journal of Physics: Conference Series (JPCS)[END_REF]12].

Actinide collisions have been studied with the TDHF approach [67,11,12,68]. Fig. 12 shows density evolutions in 238 U+ 238 U central collisions at E c.m. = 900 MeV for different initial conditions. We see that the initial orientation of the nuclei plays a crucial role on the reaction mechanism [11], with the production of a third fragment in the tip-tip collision (left), or net mass transfer in the tip-side configuration (middle). In the latter case, ∼ 6 protons and ∼ 11 neutrons, in average, are transferred from the right to the left nucleus, corresponding to the formation of a 255 Cf primary fragment.

Similar calculations have been performed on the 232 Th+ 250 Cf system [12]. An example of density evolution is shown in the right panel of Fig. 13. In this case, we observe a net transfer of nucleons from the tip of the 232 Th to the side of the 250 Cf, corresponding to an inverse quasi-fission process, i.e., the exit channel is more mass-asymmetric than the entrance channel. Indeed, in this case, a 265 Lr fragment is formed in the exit channel. It is worth mentioning that these calculations predict inverse quasi-fission for this specific orientation only, i.e., when the tip of the lighter actinide is in contact with the side of the heavier one. Indeed, the other orientations induce "standard" quasi-fission [12]. Note also that another inverse quasi-fission mechanism is predicted in actinide collisions due to shell effects in the 208 Pb region [69,6,[START_REF] Zagrebaev | Journal of Physics: Conference Series (JPCS)[END_REF]].

In the previous example, the 265 Lr heavy fragment indicates the average N and Z of a distribution. The fluctuations and correlations of these distributions have been computed with the Balian-Vénéroni prescription [32] using the tdhf3d code [70,68]. Fig. 13(left) shows the resulting probabilities assuming Gaussian distributions of the form where σ N,Z are the standard deviations and 0 ≤ ρ < 1 quantifies the correlations between the N and Z distributions [START_REF] Simenel | [END_REF]. We see that many new β-stable and neutron-rich heavy nuclei could be produced if their excitation energy is low enough to allow their survival against fission. This inverse quasi-fission process needs further investigations, in particular to predict realistic production cross-sections.

P (z, n) = 2πσ N σ Z 1 -ρ 2 -1 exp - 1 1 -ρ 2 n 2 σ 2 N + z 2 σ 2 Z - 2ρnz σ N σ Z , (12)

Conclusions and perspectives

The TDHF theory has been applied to the study of heavy ion collisions at energies around the Coulomb barrier. Its time-dependent nature allows to investigate dynamical effects responsible for the modification of the fusion thresholds. In particular, the coupling between the relative motion and proton transfer in 16 O+ 208 decreases the barrier by reducing the Coulomb repulsion. Particle transfer probabilities are predicted using a particle-number projection technique. A comparison with experimental data shows the importance of pairing correlations on transfer. The latter could be studied with new time-dependent Hartree-Fock-Bogoliubov codes. In particular, one could answer the question on the origin of these pairing correlations: Are they present in the ground-states or are they generated dynamically during the collision?

For heavy and (quasi-)symmetric systems, a fusion hindrance is observed due to the quasifission process. The possibility to study the quasi-fission mechanism with a fully microscopic quantum approach such as the TDHF theory is promising. It will help to understand the strong fusion hindrance in quasi-symmetric heavy systems and may provide a guidance for new fusion experiments with exotic beams.

Actinide collisions have been investigated both within the TDHF approach and with the Ballian-Vénéroni prescription. A new inverse quasi-fission mechanism associated to specific orientations was found. This mechanism might help to produce new neutron-rich heavy nuclei. A systematic investigation of this effect is mandatory to help the design of future experimental equipments dedicated to the study of fragments produced in actinide collisions.

Figure 1 .

 1 Figure1. Nuclear (dashed blue line) and Coulomb (dotted red line) contributions of the nucleus-nucleus potential (solid line) obtained with the HF-frozen approximation in the16 O+208 Pb system. The Wong potential[43] is shown in green dot-dashed line. The latter is obtained with a potential depth V 0 = 70 MeV, a potential diffuseness a = 0.48 fm, and nuclear radii R i = 1.25A 1/3 i fm.

Figure 2 .

 2 Figure 2. Experimental fusion barrier distribution of the 16 O+ 208 Pb system from Ref. [47]. The HF-frozen, Bass and TDHF barriers are indicated by arrows.

Figure 3 .

 3 Figure 3. Relative distance between fragments as a function of time for head-on 16 O+ 208 Pb collisions.

Figure 4 .

 4 Figure 4. (top) Density evolution for the reaction 16 O+ 208 Pb corresponding to a head-on collision at a center of mass energy E c.m. = 74.44 MeV (just below the fusion barrier). The red surfaces correspond to an iso-density at half the saturation density (ρ 0 /2 = 0.08 fm -3 ). Each figure is separated by a time step of 135 fm/c. Time runs from left to right. (bottom) Same at E c.m. = 74.45 MeV, i.e., just above the fusion threshold.

Figure 5 .

 5 Figure 5. Neutron (circles) and proton (squares) number probability distributions of the lightest fragment in the exit channel of a head-on 16 O+ 208 Pb collision at E c.m. = 74.44 MeV (left) and 65 MeV (right). Adapted from Ref. [54].

Figure 6 .

 6 Figure 6. Proton number probability as function of the distance of closest approach in the small outgoing fragment of the 16 O+ 208 Pb reaction. TDHF results are shown with lines. Experimental data (open symbols) are taken from Ref. [50].

Figure 7 .

 7 Figure 7. Distance between the centers of mass of the fragments as a function of time in head-on 90 Zr+ 124 Sn collisions for different center of mass energies.

Figure 8

 8 Figure 8. TDHF fusion thresholds for several heavy systems are compared with the proximity barrier[60] and with results from the extra-push model[START_REF] Swiatecki | [END_REF].

Figure 9 .

 9 Figure 9. Density profile in the 90 Zr+ 124 Sn head-on collision at E c.m. = 235 MeV. From Ref. [61].

Figure 10 .

 10 Figure 10. Snapshots of the TDHF isodensity at half the saturation density in the 40 Ca+ 238 U system for different initial orientations and E c.m. .

Figure 11 .

 11 Figure 11. TDHF calculations of the mass of the fragments (top) and of the quasi-fission time (bottom) in40 Ca+ 238 U central collisions as a function of the center of mass energy (divided by the proximity barrier[60]). For quasi-fission times larger than 23 zs, only a lower limit is given. Two different orientations of the 238 U are considered (see inset).

Figure 12 .

 12 Figure 12. Snapshots of the isodensity at half the saturation density in 238 U+ 238 U central collisions at E c.m. = 900 MeV from TDHF calculations. Snapshots are given at times t = 0, 1.5, 2.7, and 4.2 zs from top to bottom. From Ref. [11].

Figure 13 .

 13 Figure 13. (right) Snapshots of the isodensity at half the saturation density in 232 Th+ 250 Cf central collisions at E c.m. = 916 MeV. (left) Gaussian distributions of N and Z with widths and correlations computed with the BV prescription (linear color scale). The solid line represents the predicted β-stability line.

Acknowledgements D. Hinde, M. Dasgupta and M. Evers are warmly thanked for useful discussions. A. W. is grateful to CEA/Saclay, IRFU/SPhN, where part of this work has been performed. The TDHF calculations were performed on the NCI National Facility in Canberra, supported by the Commonwealth Government. Support from ARC Discovery grants DP06644077 and DP110102858 is acknowledged.