Wilfried Despagne
email: wilfried.despagne@icam.fr

Emmanuel Frénod
email: emmanuel.frenod@univ-ubs.fr.

MathematicS In Action Version de travail Transport hub flow modelling

Keywords: Mathematical Modelling, Ordinary Differential Equations, Haulage Activity, Hub, Express, Supply Chain Modelling. Math. classification: 00X99

The purpose of this paper is to investigate the road freight haulage activity. Using the physical and data flow information from a freight forwarder, we intend to model the flow of inbound and outbound goods in a freight transport hub. Approach: This paper presents the operation of a road haulage group. To deliver goods within two days to any location in France, a haulage contractor needs to be part of a network. This network handles the processing of both physical goods and data. We will also explore the ways in which goods and data flows are connected. We then build a first model based on Ordinary Differential Equations which decrypt the flow of goods inside the hub and which is consistent with available data. This first model is designed to work at a fine-scale level. A second model, which aggregates factors of the fine-scale model, is also built. Tests are carried out to show the accuracy of the models. Finally, an explanation on how to use the models for industrial process optimizing is given.

Introduction

This paper is part of a research program whose target is to build a model describing the workings of a fresh and frozen product transportation network.This program is funded in part by the French group STEF-TFE 1 . This group is the European leader in temperature-controlled-logistics solutions. In 2009, the company's European logistic network included 215 temperature-controlled sites and a fleet of 2,500 vehicles, transporting 30,000 tons of goods daily on over 80,000 waybills throughout Europe. In 2010, the transport branch of the group began a rigorous development process of a hub simulation system. The goal of such a system is to visualize how goods are moving in a hub under several business parameters. The target is to make it easier to fix some of those parameters used for controlling complex warehouses. Using the resulting system, we can simulate alternative cross-docking user stories, workforce inventory, goods dispatching rules and transport plans.

The nodes of the supply chain network are the hubs where merchandise is brought in, offloaded, handled, and loaded in trucks. These last 3 actions make up the cross-docking. The role of cross-docking is to move goods just-in-time, from the receiving dock to the shipping dock without creating stock, [START_REF] Udauy | Effective Cross Docking for improving Distribution Efficiencies[END_REF]. Goods are regrouped, on the hub, consignments are shipped by several senders and are intended for several customers. The hubs are connected by the comings and goings of trucks. Since the goal is to reduce the resources needed to carry goods, the transport domain is being organized into freight forwarder networks. These networks consist in several haulage contractors pooling their resources and working together to satisfy the customers' requirements. The strength of a freight forwarder network resides in the players' ability to coordinate their strengths in the user's service and then realizing the "Strategic Management of Supply and Logistics", [START_REF] Beaulieu | Gestion stratégique de l'approvisionnement et de la logistique[END_REF]. The creation of freight forwarding networks makes it possible to offer users nationwide delivery in less than 48 hours. Such a network requires strict piloting through a set of procedures and software for the management of data flow, physical streams and interfaces between the various players. This set of tools is known as Supply Chain Management (SCM). To achieve control of the flows (both goods and data), the process must be industrialized: the transport plans optimized and the flows placed in synergy through computerized data processing. The SCM tools have considerably helped carriers to improve their productivity. The carrier's key tool is TMS (Transport Management System) [START_REF] Helo | Logistics information systems: An analysis of software solutions for supply chain co-ordination[END_REF]. It makes it possible to schedule and optimize routes, carry out day-to-day operations, dialogue with users and, through the reporting function, to analyze the carrier's capacities. To be manageable, the physical flow of goods must be accompanied by a flow of data. Every package, every transfer or handling operation is recorded in databases. TMS uses the sum of this data to make information both reliable and relevant for operation and decision management [START_REF] Picot | The Fading Boundaries of the Firm: The Role of Information and Communication Technologies[END_REF][START_REF] Stefansson | Tracking and tracing: principles and practice[END_REF].

Research has been done to reduce transport costs. Efforts have brought out algorithm optimizing transport plans (see [START_REF] Bowersox | Logistical managements: the integrated supply chain process[END_REF][START_REF] Patterson | Diffusion of Supply Chain Technologies[END_REF][START_REF] Scott | Integrating the warehousing and transportation functions of the supply chain[END_REF][START_REF] Kirschfink | Intelligent Traffic Management Models[END_REF]) and based on these algorithms, software has been developed (see [START_REF] Yang | Simulation Laboratory for Evaluating Dynamic Traffic Management Systems[END_REF][START_REF] Helo | A web-based logistics management system for agile supply demand network design[END_REF]). Methods and algorithms have been set out to model and to optimize the route of goods within the hub, using operational research tools like graph theory and high dimension optimization algorithms (see [START_REF] Rohrer | Simulation and cross docking[END_REF][START_REF] Bartholdi | Reducing labor costs in an LTL crossdocking terminal[END_REF][START_REF] Iris | Positioning of goods in a cross-docking environment[END_REF][START_REF] O'kelly | Routing Traffic at Hub Facilities[END_REF][START_REF] Yu | Scheduling of inbound and outbound trucks in cross docking systems with temporary storage[END_REF][START_REF] Liu | Simulation-based personnel planning for materials handling at a cross-docking center under retail distribution environment[END_REF]). Yet, no system offers the possibility to simulate hub working in usual layout or modifications of hub characteristics or rare events impacting the hub, incorporating factors concerning workforce and material used. No system offers the possibility to optimize hub working, incorporating those factors either. In the same way, very poor possibilities are offered to simulate or optimize how hubs interact with the network they are part of. These seem to be the two topics in which improvement can be made to bring about new productivity gains.

The goal of this paper is to contribute to the understanding of the internal working of hubs, building a model and using it. We model it at a fine-scale and at a larger scale. The paper is organized as follows. In section 2, we build a hub model using ODE (Ordinary Differential Equations) type tools. We first formalize hub workings. Secondly, we build a finescale model and then, aggregating parts of it, we obtain a model, which is valid at a larger scale. In section 3, we build an implementation of the fine-scale models and we simulate the working of a simplified hub. This shows that the model fits hub qualitative functioning well. Then, the aggregated model is also implemented. Tests show a good accuracy with fine scale results. Finally, we simulate an almost realistic hub using the fine-scale model. In the "conclusion and perspective" section, we give a list of projects we plan to reach. In particular, we develop the ones concerning the way to use the fine-scale and aggregated models for optimization of business process purposes and for the building of a hub network.

Modelling of the flows transiting through the hub

This section is devoted to model building.

The Hub

We begin by formalizing hub and hub network workings. Let us look at an example used by [START_REF] Beaulieu | Gestion stratégique de l'approvisionnement et de la logistique[END_REF]. According to Fig. 1, suppliers F1, F2, etc. can each provide a different product. At the other end of the line, customers C1, C2, etc. can all want various quantities of each product. For instance, customer C1 can want 2 units of product F1, 3 units of F2, 5 units of F4, while customer C2 can want 5 units of F1, 2 units of F2 1 unit of F4 and 3 units of F5. The use of the hub enables the contractor to reduce the number of trips from 5 × 5 to 5 + 5 (Fig. 1). Besides, thanks to the hub, the contractor can offer new services, such as the processing of orders from various suppliers. Though it is not the hub's primary vocation, it can be used as a temporary storage area to offer customers extra services, such as more fluid shipping and distribution. A hub's size varies with the importance of the flow to be accommodated. The overall surface, Figure 1. The economical effect of a middleman [START_REF] Beaulieu | Gestion stratégique de l'approvisionnement et de la logistique[END_REF]. including hub, office space and access, and parking area covers an average of 1 ha. On both sides of the hub are a series of loading bays, used to load and unload trailers. A good organization of the hub is of primary importance to optimize the flows. The greatest challenges are managing available space, deadlines, quality of delivery, and safety. The hub includes several areas (Fig. 2). The "Load" and "Unload" areas face each other along the outer walls. Between them is a traffic area in which the packages are stacked on pallets and moved with forklifts and various handling apparatus. When a truck arrives, the driver hands his paperwork in to the operations department. This service assigns him an unloading bay on the "In Area" side. If the bay is already in use, the driver parks in a waiting zone. The data capturing of the transport documents is transmitted to the hub, generating the printing of labels. These labels show the consignment's origin, destination, quantity, nature, and especially the entry and exit bays it will use -and, if relevant, the "transit area" in which it is to be placed. This means that upon receipt of a consignment, an employee checks the provenance and quantity of the unloaded goods and labels them. Referring to these labels, the handlers know what route the consignment must take through the hub, and its movements are limited and controlled. In the background, the TMS regroups the packages according to their destination. If the packages do not require any sorting, they are directly taken to an area facing their loading bay. This is cross-docking. This process increases the available storage space and makes it possible, at any time, to know where the goods are, thus limiting losses. If the unloaded goods need to be sorted before departure, they are taken to the sorting zone or "Residence Area". The dispatching operations must not slow down the traffic of other packages.

To conform to the customer's required deadline, the absolute must is to have the truck depart on time. The schedule defined by operations assigns arrival and departure times for each route. If these times are disregarded, the whole day's organization suffers. This is why the trucks leave on time even when they are not fully loaded. When everything is going as planned, the loading or unloading of a truck takes about 20 minutes. The hub's productivity is calculated on a ton-per-hour index. This index, with the addition of business forecasts in tons and number of waybills, makes it possible to schedule the hub's human resources for each day of the week, [START_REF] Despagne | A Forecasting Support System for Temperature-Controlled Transport[END_REF][START_REF] Despagne | A Forecasting System Developed under R, Dedicated to Temperature-Controlled Goods Hauling[END_REF].

The position is a business vocabulary to describe the digital information to track the goods from its sender to its final destination. One position matches with one waybill. It displays data on:

• the prime contractor, who pays for transportation,

• the sender, who ships the consignment,

• the loading point: the actual place where the goods are to be collected, for example the manufacturer's warehouse,

• the final receiver of the goods, for example the retailer,

• the unloading point: the actual place where the goods are to be unloaded, for example a retailer's dispatching base,

• the refrigeration status: fresh foods, frozen foods, non-refrigerated goods,

• the nature of the goods, for example cold meats or dairy products, etc.

• miscellaneous data: delivery constraints, description of packaging, etc.

This description done, a model of hub functioning taking the above explanation into account must be set out.

Modelling procedure

We chose to follow a modelling protocol which yields models involving ODE (Ordinary Differential Equations) for several reasons. The first one is the ability of this kind of tool to model transfer of quantity with a rigorous mathematical approach. The second one is that the resulting models are flexible and well adapted to data at disposal. The third one is the fact that ODE framework offers a wide range of efficient tools for mathematical analysis. We do not tackle this topic in the present paper, but we plan to do it in forthcoming ones. The fourth one is that the resulting models can be tested via ODE based software programs. The fifth one is that it generates continuous-in-time models that can be discretized, using numerical analysis methods, at any time step, if discrete-in-time models are needed. Here is a nice advantage of our approach over direct discrete-in-time modelling: once a discrete-in-time model is gotten from the continuousin-time model, if another instantiation of the discrete-in-time model is needed with another time step, the modelling protocol does not have to be done again, only the discretization. The last reason we evoke is that optimization tools based on ODE are numerous and efficient; and we plan to use the models to optimize hub working in future works. In particular, we plan to tackle questions linked to workforce optimization, workforce allocation optimization and hub operating process optimization.

In a first step, we consider that the hub may be well described by a multi-compartment model between which mass of goods is transferred. Considering a given position (see definition on page 4), we assume that the transfer of mass of goods of this position is well quantified by transfer functions involved in an Ordinary Differential Equation system [START_REF] O'regan | Existence Theory for Non-Linear Ordinary Differential Equation[END_REF][START_REF] Arnold | Equations différentielles ordinaires[END_REF]. This system makes up what we call the fine-scale model working at the position scale. Then going over the position quantities involved in the fine-scale model, and making assumptions to deduce average transfer functions, we set out what we call the aggregated model.

Fine-scale model

We model hub activity by a five compartment model and four transfer functions. Schematically, a hub is made of five areas (Fig. 3). The first one, we call "In Area", or "In" in short, is not physically located. It takes into account that when a lorry arrives at a given hub, several tasks need to be done between registration and unloading (platform reaching, trailer opening, unloading). Goods are unloaded on an "Unload Area" or "Ulo". After unloading, goods are handled. They may be taken into a "Reserve" or to a "Storage Area" if they are scheduled to remain at the hub for a long period of time. We call the set of both areas "Residence Area" or "Res". After residing in the "Residence Area", goods are brought to a "Load Area" or "Loa" in short, before going outside which is modelled by an "Out Area" or "Out" in short. The velocity at which goods are transferred between areas depends on good packaging and on the workforce allocated to this task. Moreover, the following important facts also need to be taken into account. First, the above schematic description does not fit all the situations. For instance, goods may be directly transferred from the "Unload Area" to the "Load Area". Secondly, in big hubs, specific personnel is allocated to unloading goods, while others are allocated to handling goods within the hub, and others to loading. On the other hand, in smaller hubs, the same personnel may be allocated to any task when necessity arises.

In U lo Res Loa Out ✲ ✲ ✲ ✲ ✲ ✲ Φ In ψ In ψ U lo ψ Res ψ Loa Φ Out
Similarly, we define Φ Out,(p) (t) as the flow of mass of goods of position p leaving "Out Area" at time t. Considering quantities S A (t) = p∈P os S A,(p) (t) as the total mass of goods being in area

A at time t, we may define flow of mass between areas using a transfer function. We consider transfer function ψ In,(p) (t) which is the mass flow of goods of position p going from "In Area" to "Unload Area" at time t. We also define ψ U lo,(p) (t) as the transfer function from "Residence Area" to "Load Area" and ψ Loa,(p) (t) as the transfer function from "Load Area" to "Out Area". It is reasonable to consider that

ψ In,(p) (t) = ψIn (Π(p, t), S U lo (t), n In,(p) (t))χ(S In,(p) (t)), where χ(S) = 1 if S > 0 and 0 if S ≤ 0. (2.2) Equation (2.
2) means that ψ In,(p) (t) is a function of the position. The way that ψIn depends on p has to take into account the packaging features of the goods of position p and possibly other information concerning p. Then Π(p, t) stands for a list of information concerning p at time t. ψIn is also a function of S U lo . This allows one to take into account that S U lo may be overloaded. ψIn also depends on the number of employees n In,(p) (t) allocated to the transfer task. χ(S In,(p) (t)) says that when no more goods of position p are present in "In Area" then the transfer stops. In the same way, we also consider

ψ U lo,(p) (t) = ψUlo (Π(p, t), S Res (t), n U lo (t))χ(S U lo,(p) (t)), (2.3) ψ Res,(p) (t) = ψRes (Π(p, t), ν, S Loa (t), n Res (t))χ(S Res,(p) (t)), (2.4
)

ψ Loa,(p) (t) = ψLoa (Π(p, t), S Out (t), n Loa (t))χ(S Loa,(p) (t)). (2.5)
Having those quantities and transfer functions, we can write the ordinary differential equation system modelling the transfer of mass of goods of position p within the hub.

dS In,(p) dt = Φ In,(p) -ψ In,(p) , (2.6
)

dS U lo,(p) dt = ψ In,(p) -ψ U lo,(p) , (2.7
)

dS Res,(p) dt = ψ U lo,(p) -ψ Res,(p) , (2.8)
dS Loa,(p) dt = ψ Res,(p)ψ Loa,(p) , (2.9)

dS Out,(p) dt = ψ Loa,(p) -Φ Out,(p) , (2.10)
translating that in each area, good mass variation is the result of an in-flow and an out-flow.

Once this model is set out, we have to show that it is well adapted to the working of big hubs by fixing functions ψA and Φ A . We also have to show that it is flexible enough to take into account smaller hubs. For this, we present what we call "area overlapping" and adapted workforces n A (t).

The information package Π(p, t) is the following

Π(p, t) = (p, m(p), N p (p, t), T In (p), T Out (p), P ub(p)). (2.11)
Where m(p) is the mass of goods constituting position p and N p (p, t) its pallet number. We consider that this pallet number depends on time t to be able to take into account repackaging that may be applied to positions within the hub. T In (p) is the time at which the truck containing position p is ready to be unloaded and T Out (p) is the time at which position p has to leave the hub. P ub(p) is the list of positions which are in the same truck as p and that have to be unloaded before p. In big hubs, the following form may be chosen for ψ In :

ψIn (Π(p, t), S U lo (t), n In,(p) (t)) = 1 τ In m(p) N p (p, t) n In,(p) (t)χ(S U lo max -S U lo (t)). (2.12)
In this equation, τ In stands for the time needed to transfer a pallet from "In Area" to "Unload Area" and S U lo max for a maximal capacity. Hence, since m(p) Np(p,t) is the mean mass per pallet, expression (2.12) says that the mass per pallet from "In Area" to "Unload Area" is in direct proportion to the mass of each pallet and to the workforce allocated to unload position p, and that it is in inverse proportion to the time needed to unload one pallet. Factor χ(S U lo max -S U lo ) makes the mass flow stop when S U lo max is reached. Notice that S U lo max may be set to +∞ if this threshold effect is not necessary. Workforce n In,(p) (t) allocated to unload position p may be defined as a proportion, linked with the mass fraction of position p in the total mass to unloaded modulated by a priority index, of the total workforce n In (t) allocated to unloaded task by the following formula:

n In,(p) (t) = t t-∆t γ(p, s)S In,(p) (s)ds t t-∆t p ′ ∈P os γ(p ′ , s)S In,(p ′ ) (s)ds n In (t).
(2.13)

Taking time average in this formula permits to allocate a large workforce to a large mass, but only in mean. The time interval length ∆t, on which the values are averaged, is the relaxation time defined in page 6. In (2.13), γ(p, t) ∈ [0, +∞] is the priority index. The very target of this priority index is to enforce the workforce allocated to a position with a high value of priority index so that it is processed faster. The question of the definition of γ(p, t) is important. A simple choice consists in setting

γ(p, t) = γ(p) = γ T Out (p) -T In (p) , for any t ∈ [T In (p), T Out (p)]
and for a given constant γ,

which allocates a high value to priority index to positions that are planned to pass quickly through the hub. A second reasonable definition consists in setting a time dependent function

γ(p, t) = min(γ 1 + γ2 T Out (p) -t ; γ3 ), for any t ∈ [T In (p), T Out (p)] (2.15)
for three given constants γ1 , γ2 and γ3 . This choice causes an enforcement of the priority index when time increases.

Remark -It may be noticed that in a framework of hub operating process optimization, γ(p, t) may be considered as an unknown function that has to be set via an optimization under constraint process.

Before defining the other functions ψA for A ∈ {U lo, Res, Loa}, a notation generalization needs to be set: τ U lo stands for the time needed to transfer a pallet from "Unload Area" to "Residence Area", τ Res stands for the time needed to transfer a pallet from "Residence Area" to "Load Area" and τ Loa stands for the time needed to transfer a pallet from "Load Area" to "Out Area". S A max is the maximal capacity of area A and n A (t) and n A,(p) (t), linked by

n A,(p) (t) = t t-∆t γ(p, s)S A,(p) (s)ds t t-∆t p ′ ∈P os γ(p ′ , s)S A,(p ′ ) (s)ds n A (t), (2.16)
are the workforces allocated to area A and the workforce allocated to position p in area A.

Then, the following forms may be chosen for ψUlo and ψLoa :

ψUlo (Π(p, t), S Res (t), n U lo,(p) (t)) = 1 τ U lo m(p) N p (p, t) n U lo,(p) (t)χ(S Res max -S Res (t)), (2.17
)

ψLoa (Π(p, t), S Out (t), n Loa,(p) (t)) = 1 τ Loa m(p) N p (p, t) n Loa,(p) (t)χ(S Out max -S Out (t)). (2.18)
Function ψRes,(p) depends on a supplementary variable ν which is the time at which the transfer of goods of position p, from "Residence Area" to "Load Area", has to begin. This variable helps take into account that goods may be scheduled to remain a long period of time within a given hub. The definition of ψRes,(p) is then:

   ψRes,(p) (Π(p, t), ν(p), S Loa (t), n Res,(p) (t)) = 1 τ Res m(p) Np(p,t) n Res,(p) (t)χ(S Loa max -S Loa (t)) if t ≥ ν(p), = 0 otherwise. (2.19)
For the definition of γ(p, t) which is needed in n Res,(p) computation, if before arriving in "Residence Area", it was defined by (2.14), then it is necessary to set its value to

γ(p, t) = γ T Out (p) -ν(p) for any t ∈ [ν(p), T out (p)]. (2.20)
On the other hand, if before arriving in "Residence Area", it was defined by (2.15) this definition may be kept.

TRANSPORT HUB FLOW MODELLING

Concerning ν(p) several choices may be carried out. The simplest one consists in setting

ν(p) = T Out (p) -ν, (2.21)
for constant ν which is an estimated time for a given position to go from "Residence Area" to "Out Area". Of course, constant ν(Π(p, t)) may depend on Π(p, t) with the idea that the more pallets p has, the longer ν(Π(p, t)) has to be. ν(p) may also consist in an unknown that has to be set by an optimization process.

Considering Φ In,(p) and Φ Out,(p) , denoted by δ [t=T ] the Dirac mass at time T and by δε,T (t) its regularization with temporal expansion ε:

               δε,T (t) = 0 if t < T -ε, = t -T + ε ε 2 if T -ε ≤ t < T , = -t + T + ε ε 2 if T ≤ t < T + ε, = 0 if t ≥ T + ε, (2.22)
we may set

Φ Out,(p) = m(p)δ [t=T Out (p)] , (2.23)
which says that the total mass of position p is transferred out of "Out Area" at time T Out (p), or

Φ Out,(p) = m(p)δ ∆t 4 ,T Out (p) (t), (2.24)
which is a regularization of (2.23) with temporal expansion ε = ∆t/4, where ∆t is the relaxation time defined on page 6. Concerning Φ In,(p) , we introduce a time T (p) needed for unloading positions in P ub(p) before reaching p in the truck, and then we define

Φ In,(p) = m(p)δ [t=TIn(p)+ T (p)] , (2.25)
or

Φ In,(p) = m(p) δ ∆t min 4
,T In (p)+ T (p) (t).

(2.26)

In those definitions, T (p) may be chosen as

T (p) = τ In p ′ ∈P ub(p) N p (p ′ , T In (p)).
(2.27) It may also be set via an optimization process that we will see below.

The model we have just built is well adapted to big hubs. For smaller hubs, or for hubs which do not work exactly as described above, modifications and simplifications of the model may be done in order to adapt it. Some of those modifications will now be explained.

First, the workforce may not be allocated by area, but by area group. In this situation, A stands for the area group. For instance, A may be {In, U lo} or {Res, Loa} or {In, U lo, Res, Loa}. Then, for a given A which is in group A, the following formula Notice that it is possible to overlap more than two areas and that this "area overlapping" is compatible with any workforce distribution of type (2.28) or (2.16).

n A,(p) (t) = t t-∆t min γ(p, s)S A,(p) (s)ds t t-∆t min p ′ ∈P os A ′ ∈A γ(p ′ , s)S A ′ ,(p ′ ) (s)ds n A (t), (2

Aggregated Models

We will now build an aggregated model which is valid regarding flows more macroscopically and at a larger time scale. (Typically, the minimum time step of a discrete time model deduced from the aggregated model is about 4∆t min ≃ 1h). The aggregated model involves less variables and data and is easier to implement. Hence, having it at hand may be important, for instance, for the first operational implementations and for network modelling. The main characteristic of this model is that it is not as precise as the fine-scale one defined by equations (2.1)-(2.27). In particular, it aggregates flows from every position.

We first consider that function Φ In (t) gives the flow entering "In Area" at time t and Φ Out (t) gives the flow leaving "Out Area" at time t. Then the quantities involved are:

S A (t) = p∈P os S A,(p) (t) for A ∈ {In, U lo, Res, Loa, Out}, (2.30)
where the sum may be reduced to p ∈ P os such that T In (p) ≤ t ≤ T Out (p).

We deduce the aggregated model from equation ( ) is built to be conform to the functioning of big hubs. Nonetheless, for smaller hubs "area overlapping" in the same spirit of the process described by, (2.29) may be done.

Validation tests

The model described in section 2.3 and 2.4 was programmed in environment R [START_REF]R: A Language and Environment for Statistical Computing[END_REF]. The numerical method used for computations is the Runge-Kutta 4 methods, [START_REF] Butcher | Numerical methods for ordinary differential equations[END_REF]. The following sections shows that the models are applicable in real situations. At first, basic user stories are applied to highlight, under different parameters, the flow generated by two positions. Simulations are done with booth models. At the end of the section, a stronger user story simulates a more realistic situation with thirty waybills spread over 30 hours of cross-docking.

Fine-scale model results

To show the model's performances, we simulated the activity of a hub from a simple example. We had two positions come within the space of an hour. They stayed at the hub for ten hours. We wanted to observe how their mass filled the different hub areas during these ten hours.

To do so, the above-described model requires set parameters for positions, available workforce at the hub, and the size of the hub. Four user stories will help us understand how a model behaves according to the varied parameters. The first user story sets the parameters as indicated in the following tables. Two positions of 5 tons and 20 pallets each come into the hub, respectively at 1:00 and 2:00. They are programmed to leave at 11:00 and 12:00. Their parameters to determine the priority index with formula (2.14) are the same. The positions can be unloaded at the same time. The different hub areas are large enough to store the total mass of both positions. The amount of manpower is equally divided between each position. The transfer time of a pallet from one area to another is estimated at one minute. To estimate the position transfer time from the "Residence Area" to the "Out Area", we ran the model once to see the time necessary for unloading. This amount of time is allotted to parameter ν.

Figure 4 shows the results obtained by running the model with the above parameters. The first graph in Fig. 4 shows the evolution of the sum of mass of goods from both positions, in each area. It also shows the evolution of the total mass stored in the hub. The second and third graphs in Fig. 4 show the evolution of the mass in each area, for each position. The "In" curve represents the evolution of mass of goods in the "In Area". The "Ulo" curve represents the evolution of mass of goods in the "Unload Area". The "Res" curve represents the mass of goods in the "Residence area". The "Loa" curve represents the evolution of the mass of goods in the "Load Area". The "Out" curve represents the evolution of mass of goods in the "Out Area". While the "Unload", "Residence", and "Load" areas are physically marked areas in the hub, the "In" and "Out" areas represent the sum of the areas in the trailers which are docked at the hub for unloading and loading. Indeed, a truck contains goods from a series of positions. The truck and its trailer go to a dock to be unloaded at a specific time. After being handled and moved from the "Unload Area" to the "Residence Area", then from the "Residence Area" to the "Load Area", the goods from the same position are loaded in the trailer for departure at a scheduled time. This is the date of departure or the date at which the goods leave the "Out Area". The "In" curve from the first graph in Fig. 4 shows that the start date for the unloading of position n°1 is 1:00. The second graph shows that the unloading continues to approximately 7:30. This graph also shows that the transfer of position n°1 mass is slowed down from the "In Area" to the "Unload Area". This is simply due to the coming into the hub of the second position. The evolution of the distribution of second position (n°2) mass throughout the hub can be seen in the third graph in Fig. 4. The second position comes in for unloading at 2:00 and leaves the "Out Area" at 11:00. As previously stated, the manpower available for unloading is equal to the manpower available for transferring goods from the "Unload Area" to the "Residence Area". Thus, as soon as the goods are unloaded, they are transferred to the "Residence Area". This is why, as can be seen in the three graphs from Fig. 4, the "Unload Area" does not keep stock. While the "In Area" is being cleared, the "Residence" and "Out" areas are filling up (see "Res" and "Out" curves). The "Out Area" starts filling up at 2:00 because we set parameter ν at 8:00. The "Load Area" does not fill up. The reason for this is that the goods leaving the "Resident Area" are immediately transferred to the "Out Area". Position n°1 leaves, as scheduled, at 10:00 and position n°2 leaves at 11:00. When the goods go to the "Out Area", the "In" and "Resident" areas empty out. The "Tot" curve shows the goods mass present in the hub at each moment.

The second user story uses the same parameters as the previous one except that the available manpower to transfer the goods from "Unload Area" to "Resident Area" is reduced. The amount of manpower goes from 0.1 to 0.08. Fig. 5 gives the results obtained. As before, the "In" curve has two peaks which correspond to both positions coming into the "In Area". The drops in the curve between 1:00 and 2:00 and then between 2:00 and 7:00 correspond to the unloading times of the trucks. As the team to transfer the unloaded goods from the "Unload Area" to the "Resident Area" is smaller than the unloading team, the goods progressively accumulate in the "Unload Area" (see "Ulo" curve). This is also why the "Resident Area" does not fill up as much as in the previous situation. It can also be noted that the "Res" curve starts dropping at 3:00. This is due to the fact that the goods transferred to the "Resident Area" are transferred to the "Load Area" right away, to be moved to the "Out Area". There is a break in the "Out" curve slope at 3:00, being the time at which the second positions start being loaded into the trucks. The "Out" curves reach the level of 10 tons at 9:00. This is the sign that the goods in both positions are loaded and ready to leave the "Out Area". At 10:00, the "Out" curve drops for the first time to show the departure of position n°1. Then, the curve stagnates at 5 tons at 1 hour and drops again at 11:00 to 0 ton. Position n°2 is left.

Similarly, when the amount of available manpower decreases for loading (see parameter n Loa ), goods stock pile up in the "Load Area", therefore causing the "Loa" curve to be positive.

A third user story sets the parameters as they were for the first simulation, considering that the "Residence Area" is not large enough to handle more than one ton of goods (S Res max = 1). The results illustrated in Fig. 6 indicate that the consequence of this constraint is the stocking of surplus goods in the "Unload Area". Fig. 6 perfectly shows that the "Resident Area" (see "Res" curve) starts filling up at 1:00 with the arrival of the first position. At 1:30, it reaches its maximum set capacity of one ton. At the same time, the "Unload Area" begins to stock pile goods. The break in the "Ulo" curve at 2:00 indicates that the goods are coming from the first and second positions. As expected, the "Unload Area" empties out before the "Resident Area".

The last user story considers that six and one half hours to unload two positions of 5 tons is long. Six and one half hours is how long the "In" curve stays positive, thus between 1:00 and 7:30. The first idea that comes to mind is to decrease the unloading time and to increase the amount of manpower allotted to unloading. We set all of the parameters as they were in the first user story, but we increased the value of parameter n In from 0.1 to 0.5. The "In" curve in Fig. 7 indicates that the positions were unloaded faster than in the first user story. Both were unloaded in less than an hour. However, we did not increase the amount of manpower for transferring the goods from the "Unload Area" to the "Resident Area". This is why it can be noted that the "Ulo" curve grows, which means that the goods stay in the "Unload Area" longer. Again, the "Unload Area" becomes empty long before the "Resident Area", which indicates that the goods are constantly being moved. Furthermore, the sum of the masses in the hub is 5 tons between 1:00 and 2:00, and between 10:00 and 11:00. It is 10 tons between 2:00 and 10:00. This curve is an indicator of the validity of the model. These four user stories show that, for these simple cases, the model behaves in accordance with the express hub. They also show that the model, as it was programmed, allows one to simulate several user stories. The model shows, through irregular values (negative ones, for example) that a given user story is not possible. For example, Fig. 8 shows the result when we use the first user story and have the first position leave at 6:00. The model tells us, through negative values for the "Out Area", that this objective is not feasible with the means used. To achieve it, the capacity of manpower allotted to this position would have to be increased throughout the various hub areas. The distribution of weight throughout the hub can be monitored position by position, therefore allowing us to easily detect any problems.

Results obtained for the aggregated model

Like the fine-scale model, the aggregated model was developed in an R environment. The goal was not to program it for industrial use, but to show the credibility of the mathematical model. We would also like to show that the response times of the aggregated model are faster than those of the fine-scale model. To demonstrate the similarities of results obtained from both models, we ran the aggregated model with parameters identical to the above-described stories (section 3.1), and others which are the results of the aggregation of parameters of the fine-scale model. Finally, we superposed the curves from both models. The aggregated model requires fewer parameters than the fine-scale model. It does not give the details of the distribution of goods throughout the hub, position by position. To keep the same user story as the first one given in section 3. Figure 9 superposes the results obtained by both models studied. It shows that the models empty out the "In Area" at the same speed. The transfer of the goods from the "Unload Area" to the "Resident Area" also takes place in the same way for both models. Both models start loading at the same time. However, the fine-scale model loads the "Out Area" a little more slowly than the aggregated model. This is why the "Resident Area" does not fill up as much with the aggregated model. This phenomenon could come from function n A,(p) (t) (2.16) which does not exist for the aggregated model. As expected, the "Out Area" is cleared at 10:00 and 11:00. In reality, this is the time at which the trucks leave.

In the same conditions, the run time of the first user story with the fine-scale model was 4:47min while the run time of the aggregated model was 1:37min, three times less than the former.

A heavier test

In this section, we will show that the fine-scale model works in a realistic situation. Schematically, STEF-TFE hub network working may be described in the following way (see [START_REF] Despagne | A Forecasting Support System for Temperature-Controlled Transport[END_REF]): Day begins by "collection rounds". After a "collection round", truck brings goods to the receiving dock between 8 a.m. and noon. The goods are cross-docked through the hub according to their destinations and shipped between 6 p.m. and 10 p.m. to a "consolidation round" or between 8 a.m. and noon to a "delivery round". After "consolidation round", trucks bring goods to the receiving dock between midnight and 5 a.m. The goods are cross-docked through the hub according to destination and are shipped between 8 a.m. and 12 p.m. to a "delivery round" or between 6 p.m. and 10 p.m. to a "consolidation round". "Collection round" is the action consisting in taking a vehicle to the sender's dock to load the goods and carry them to the receiving dock. "Consolidation round" is the transport operation consisting in sending a consignment between two hubs. "Delivery round" is the transport from the dispatching hub to the final receiver.

Consider that we want to simulate the flow of goods inside the hub over a 30 hour period. This will allow us to analyze the consolidation and delivery activities. The first one is done by the day team (from about 8 a.m. to about 6 p.m.), and the second one by the night team (from about 12 p.m. to about 8 a.m. the next day). We inject for the consolidation activity 131 tons split into 20 waybills. We also inject for the delivery activity 143 tons split into 28 waybills. We fix that a pallet contains one ton. We set the max weight of each area to 100 tons and the manpower available for each area to 1.2. the other hub parameters are the same as the example of section 3.1. Figure 12 shows the position parameters used to simulate the flow of goods through the hub over a period of 30 hours. Between 6 p.m. and 12 p.m. of the first day, the activity is low. Trucks are leaving hub. On the contrary, between 6 a.m. and 8 a.m. trucks arrive at the receiving dock every half hour. At 6 a.m. warehousemen have to unload 28 tons (5 + 8 + 2 + 3 + 10), at 6:30 a.m. they have to unload 28 tons (5 + 15 + 2 + 3 + 3) and so on. This is well illustrated in figure 12. The goods are loaded for consolidation rounds that have to leave the hub between 6 p.m. and 10 p.m. Every hour, trucks leave the shipping dock. The first departure is at 6 p.m. with 31 tons (3 + 10 + 15 + 3) of merchandise. Figure 10 shows the results of the simulation. We perfectly note the two activity levels described above. Between 6 a.m. and 8 a.m of the first day, trucks arrive from collection rounds to the receiving docks. Trucks are quickly unloaded. Goods are transferred in the "Residence Area". They stay there between 8 a.m. and noon. From 6 a.m. to 8 a.m., the workforce is able to transfer goods from the "In Area" to the "Residence Area" without staying in the "Unload Area". From 8 a.m. to 12:30 p.m., the "Residence Area" is full so that the "Unload Area" becomes stock. Between 11:30 a.m. and 2 p.m., goods are transferred from the "Residence Area" to the "Out Area" without staying in the "Load Area". From 2 p.m. to 7 p.m., the "Out Area" is full so that goods are stocked in the "Load Area". From 7 p.m., we observe consolidation rounds departure consequences. Between 10 p.m. and midnight, the hub is nearly empty. At midnight begins the consolidation rounds arrivals. The "In Area" receives goods every 30 minutes. These goods are transferred to "Residence Area" without staying in the "Unload Area". We can see that the "Out Area" receives goods from the truck unloading that begins. Goods are then moved from the "In Area" to the "Out Area" almost without staying in the "Unload", "Residence" or "Load Areas". Between midnight and 4 a.m. of the next day, the workforce must treat goods on the five areas, that's why a little stock accumulates on the "Residence Area". From 3 a.m., the "Out Area" is full and stock accumulates on the "Load Area". From 8 a.m., trucks leave the shipping dock for a delivery round. All those described facts are well in accordance with a real hub working. Yet, we can point two inconsistencies that will need to be remedied in future works by adding constraints in the model: The "Out Area" corresponds, in the real world, to an area located in the dock. Goods can be transferred in the "Out Area" if trucks are present and goods leave "Out Area" as the consequence of departure of some trucks. In real situations, trucks don't stay at the shipping dock five hours before they leave, as figure 10 seems to show. On another hand, at 8 a.m the second day, the hub cannot ship at once to 36 tons (2 + 5 + 5 + 1 + 5 + 7 + 5 + +2 + 4). That is observed on the curve representing the total mass in the hub. Normally, the curve should slump from 143 tons to 107 tons at 8 a.m. of the second day. Here, the curve slumps from 143 tons to 114 at 8 a.m. and takes half an hour for a decrease of 7 tons more. We can presume that the workforce is not numerous enough to deal with this task.

As we wrote above, the model was programmed with software R. However, R cannot handle a large number of positions. The above examples show the management of some positions in the hub. To manage 2 positions, our laptop, having a 1.8GHz dual processor and 2.5 Go of RAM required a little less than 5 minutes. Given that a transport agency manages between 1, 000 and 2, 000 positions a day, our computer, which is certainly not powerful, would need more than three and a half days to simulate a real life user story. We believe that, to significantly increase execution performances, one solution would be to program the mathematical model in a compiled language such as C++.

Perspectives and conclusion

Perspectives

Beyond the simulations we did to validate our models, a first perspective will be to develop a software program that can be parametrized to fit any real hub. This program will be used to simulate hubs in their standard workings to make users understand those workings. It will also be used to simulate effect of modifications of hub characteristics or of dysfunctions. Yet, possessing such a tool can also bring resolution of hub operating processes via optimization using the program. This is a long term perspective, but in order to illustrate the way we plan to explore, we give examples of how the just built models may be used in order to translate industrial problems linked with the working of the hub into optimization questions.

The first example consists in setting the question of determining the necessary workforce and its division between areas in terms of a constrained minimization problem involving the aggregated model. We consider the aggregated model (2.30)-(2.45) with known "In" and "Out" flows Φ In (t) and Φ Out (t). Inner flows ψ In , ψ Res and ψ Loa are given by (2.33)-(2.40) considering that workforces n In (t), n U lo (t), n Res (t) and n Loa (t) are unknown to be determined. At a given time t i , we consider that S In (t i ), S Res (t i ), S Loa (t i ) and S Out (t i ) are known and, a final time t f being given, we want to find functions n, n In , n Res , n Loa and n Out defined on [t i , t f ], which are linked with: The second example concerns the allocation of the workforce to positions which will be seen as an optimization process on the priority index γ(p, t) involved in equation (2.14). Here, we consider the fine-scale model (2.1)-(2.27) for every position p that passes by the hub between an initial time t i , where every S A,(p) (t i ) are known for every concerned position and area A ∈ {In, U lo, Res, Loa, Out} and a final time t f . We consider that n In (t), n U lo (t), n Res (t), n Loa (t) are given for any time t ∈ [t i , t f ]. For instance, they may have been computed by the procedure just described. We also consider that Φ Out,(p) (t) are known for each of the considered positions, for instance given by (2.23) or (2.24). Concerning Φ In,(p) (t), we consider that it is given by (2.25) in which T (p) is an unknown to be determined. Concerning the inner flows, we consider that they are given by formula (2.2)-(2.5),(2.12) and (2.17)- (2.19) where ν(p) = ν(Π(p, t)) are given and n A,(p) given by (2.16), for A ∈ {In, U lo, Res, Loa, Out} with γ(p, t) being unknown functions to determine. As a matter of fact, we are looking for constants T (p) and function γ(p, t) such that p/S A,(p) (t)>0 n A,(p) (t) ≤ n A (t), for A ∈ {In, U lo, Res, Loa, Out},

n In (t) + n U lo (t) + n Res (t) + n Loa (t) = n(t), ∀t ∈ [t i , t f ] , ( 4
γ(p, t) = 0 for all t ∈ T In (p), T In (p) + T (p) ,

where,

T (p) = τ In p ′ ∈P ub(p) N p (p ′ , T In (p)) n In,(p ′ ) (t) , (4.6)
and The third example consists in looking at the setting of times ν(p) as a constrained minimization problem. The value of ν(p) gives an estimation of the time needed for position p to go from "Residence Area" to "Out Area". We assume that the number of pallets N p (p, t) of position p remains the same from ν(p) to T out (p) and we consider that ν(p) is given (2.21) with ν(p) given by:

dS Out,(p) dt (t) = ψ Loa,(p) (t) -Φ Out,(p) (t) ≥ 0 for any, t ∈ [t i , t f ] , (4.7
ν(p) = ν(Π(p, t)) = aN p (p, T out (p)) + b, (4.9
dt (t) = ψ Loa,(p) (t) -Φ Out,(p) (t) ≥ 0 for any t ∈ [t i , t f ] . (4.11)
In the near future, we also plan to use the hub models built here to model a hub network working. We briefly explain now the way to do it. The main issue is to build connexions between instantiations of model ( 2 2.45) kind cannot be accessed. For instance, H unm may contain hubs of subcontractors. Then, we consider the sets K i,j , for i = 1, . . . , n and j = 1, . . . , n including all the possible routes from hub #i to hub #j. Notice that if hub #i is not connected with hub #j we have K i,j = ∅. It also seems reasonable to write K i,j in the following way:

K i,j = {k 1 i,j , k 2 i,j , . . . , k l i,j , . . . , k L i,j i,j }. (4.12)
Where k l i,j stands for a possible route from hub #i to #j. To each route k l i,j is attributed a time T l i,j which is the time required to travel k l i,j . Concerning the mass of goods flow of position p entering "In Area" of hub #i, Φ In,(p) i (see formula 2.1), we introduce, R λ,(p) i (t) which is the mass flow of good of position p coming from factory or depot λ, B λ,(p) i (t) which is the mass flow directly brought by factory λ at hub #i, and F j,(p) i (t) which is the mass flow of goods of position p issued from hub #j. Then we have:

Φ In,(p) i (t) = λ∈Λ i R λ,(p) i (t) + λ∈Λ i B λ,(p) i (t) + n j=1 F j,(p) i (t), (4.13)
where Λ i is the set of factories working with hub #i.

Conclusion

In this paper, we respond to a request expressed by our partner specialized in logistics. This partner wants to build a model describing the workings of a fresh and frozen product transportation network.

The first part of the research program was to predict daily inflows and outflows of goods in and from a hub [START_REF] Despagne | Construction, analysis and implentation of a forecasting model[END_REF].

The second part of this is reported in the present paper. In it, we begun by introducing how a national network of logistics hubs works. We explained what cross-docking is. We produced Ordinary Differential Equation models describing hub working at the fine scale and at a coarser one. Those models involve a large number of parameters, but protocols to set some of them are explained. So that, actually, resulting models take into account a variable number of parameters. We showed four different user stories that help to provide answers to business linked questions. We showed how express other business linked questions in optimization terms involving the models.

The third part of the research program, which needs to be led from now is to simulate how hubs interact together. Ideas to reach that goal are briefly given in the end of the present paper.

Figure 2 .

 2 Figure 2. Flow of goods within a hub. "Cross-docking" on the left and with "Residence Area" on the right.

Figure 3 .

 3 Figure 3. The five compartments and their transfer functions.

Figure 4 .

 4 Figure 4. Goods mass distribution, coming from 2 positions, in the various hub areas, with n Res = 0.1, followed by details for each position.

Figure 5 .

 5 Figure 5. Goods mass distribution, coming from 2 positions, in the various hub areas, with n Res = 0.08.

Figure 6 .

 6 Figure 6. Goods mass distribution, coming from 2 positions, in the various hub areas, with S Res max = 1, followed by details for each position.

Figure 7 .

 7 Figure 7. Goods mass distribution, coming from 2 positions, in the various hub areas, with n In = 0.5.

Figure 8 .

 8 Figure 8. Incoherence in goods mass distribution in the various hub areas.

  The model we create takes all of the previous facts into account. Moreover, it has to be valid at a time scale of several minutes. Indeed, the applications for which we want to use it need models responding with a relaxation time ∆t ranging from a minimum relaxation time ∆t min which is about 15 min to a maximum relaxation time ∆t max which is about 3 hours. Hence, we build an ordinary differential equation model. If needed for applications, from this ordinary differential equation model, discrete models may be deduced by approximation methods. For A ∈ {In, U lo, Res, Loa, Out}, p being in the set P os of positions and t ∈ R, we define S A,(p) (t) as the mass of goods belonging to position p at time t in area A. Strictly speaking, taking into account the time scale considerations S A,(p) (t) is the mean value of this mass of goods over an interval with length the relaxation time ∆t and centred in t. We also define Φ In,(p) (t) which is the mass flow of goods of position p entering "In Area" at time t. Its precise definition consists in saying that, once the relaxation time ∆t, ∆t min ≤ ∆t ≤ ∆t max is defined, the mass entering "In Area" between t and t + ∆t is

	t+∆t	Φ In,(p) (s)ds.
	t

  Res (t), which is the mass flow of goods from "Residence Area" to "Load Area", we consider that the fact that goods may remain a long time in "Residence Area" can be described by a known function S Res Rem which gives, at any time, the mass of goods which is scheduled to remain in "Residence Area". Then setting ψRes (S Loa , n Res ) = ω Res n Res χ(S Loa max -S Loa ),

	1 τ A Those factors may be averaged from history or using data or running simulation using the fine-m(p) N p (p, t) for A{In, U lo, Res, Loa, Out}. (2.31) scale model in order to deduce averaged coefficients ω A for A{In, U lo, Res, Loa, Out}. (2.32) Then, we define transfer function ψ In (t) which is the mass flow of goods going from "In Area" to "Unload Area" at time t by ψ In (t) = ψIn (S U lo (t), n In (t))χ(S In (t)), (2.33) with ψIn defined by ψ Res (t) is defined by ψ Res (t) = ψRes (S Loa (t), n Res (t))χ(S Res (t) -S Res Rem (t)). (2.40) The ordinary differential equation system finally reads dS In dt = Φ In -ψ In , (2.41) dS U lo dt = ψ In -ψ U lo , (2.42) dS Res dt = ψ U lo -ψ Res , (2.43) dS Loa dt = ψ Res -ψ Loa , (2.44) dS Out dt = ψ Loa -Φ Out . (2.45) ψIn (S (2.39) Model (2.30)-(2.45

2.1)-(2.27). In equation (2.12),(2.17),(2.18), and (2.19) we find factors: U lo , n In ) = ω In n In χ(S U lo max -S U lo ). (2.34) In a similar way, we define ψ U lo (t) = ψUlo (S Res (t), n U lo (t))χ(S U lo (t)), (2.35) ψ Loa (t) = ψLoa (S Out (t), n Loa (t))χ(S Loa (t)), (2.36) where ψUlo (S Res , n U lo ) = ω U lo n U lo χ(S Res max -S Res ), (2.37) ψLoa (S Out , n Loa ) = ω Loa n Loa χ(S Out max -S Out ). (2.38) Concerning transfer function ψ

Table 1 .

 1 Position parameters.

	Position number:	Mass	in	Number of	Positions	Unloading	Outbound
	p		tons: m(p)	pallets: N p	ahead:	date: T In (p)	date:
						P ub(p)		T Out (p)
	1		5		20	N ull	1	10
	2		5		20	N ull	2	11
	Position p transfer time	Parameter 1 of	Parameter 2 of	Parameter 3 of
	between	"Residence	the priority index:	the priority index:	the priority index:
	Area" and "Out Area":	γ 1		γ 2		γ 3
	ν
	8			1		1		5
	8			1		1		5

Table 2 .

 2 Position parameters, continued.

	Max. weight in	Max. weight in	Max. weight	Max. weight in	Max. weight in
	"In Area", in	"Unload Area",	in "Residence	"Load Area", in	"Out Area", in
	tons: S In max	in tons: S U lo max	Area", in tons: max S Res	tons: S Loa max	tons: S Out max
	20	20	20	20	20

Table 3 .

 3 Hub parameters.

	Transfer time for	Transfer time for	Transfer time for	Transfer time for
	a pallet from the	a package from	a package from	a package from
	"In Area" to the	the "Unload Area"	the	"Residence	the "Load Area"
	"Unload Area" in	to the "Residence	Area"	to	the	to the "Out Area"
	fractions of hours:	Area" in fraction	"Load Area" in	in fraction of
	τ In	of hours: τ U lo	fraction of hours:	hours: τ Loa
			τ Res
	1/60	1/60	1/60			1/60

Table 4 .

 4 Hub team parameters.

Table 5 .

 5 Hub team parameters, continued.

	Manpower avail-	Manpower avail-	Manpower avail-	Manpower avail-
	able for unload-	able for transfer	able for transfer	able for loading:
	ing: n In	to the "Residence	to	the	"Load	n Loa
		Area": n U lo	Area": n Res
	0.1	0.1	0.1			0.1
		∆t min
		1/4

Table 6 .

 6 Model parameters.

Table 7 .

 7 1, the aggregated model requires the following parameters: Position parameters.

	Position Number Mass in tons:	Unloading	Outbound	Position	p
		m(p)	date: T In (p)	date: T Out (p)	transfer time
					between "Res-
					idence Area"
					and	"Out
					Area": ν
	1	5	1	10	8
	2	5	2	11	8
	Max. weight in	Max. weight in	Max. weight	Max. weight in	Max. weight in
	"In Area", in	"Unload Area",	in "Residence	"Load Area", in	"Out Area", in
	tons: S In max	in tons: S U lo max	Area", in tons: max S Res	tons: S Loa max	tons: S Out max
	20	20	20	20	20

Table 8 .

 8 Hub parameters.

	Mean	hourly	Mean	hourly	Mean hourly hu-	Mean hourly hu-
	human cost to	human cost to	man cost to trans-	man cost to trans-
	transfer a unity of	transfer a unity of	fer a unity of mass	fer a unity of
	mass between the	mass between the	between the "Res-	mass between the
	"In Area" and the	"Unload	Area"	idence Area" and	"Load Area" and
	"Unload	Area":	and the "Res-	the "Load Area":	the "Out Area":
	ω In		idence	Area":	ω Res	ω Out
			ω U lo
	15		15		15	15

Table 9 .

 9 Hub parameters, continued.

	Manpower avail-	Manpower avail-	Manpower avail-	Manpower avail-
	able for unload-	able for transfer	able for transfer	able for loading:
	ing: n In	to the "Residence	to	the	"Load	n Loa
		Area": n U lo	Area": n Res
	0.1	0.1	0.1			0.1

Table 10 .

 10 Hub team parameters.

	∆t min
	1/4

Table 11 .

 11 Model parameters.

  .[START_REF] Udauy | Effective Cross Docking for improving Distribution Efficiencies[END_REF]

	Wilfried Despagne & Emmanuel Frénod
	and which minimizes
		t f
		n(t)ds,	(4.2)
		t i
	under the constraint
	dS Out dt	= ψ Loa -ψ Out ≥ 0.	(4.3)

  )for every position being in the hub between initial time t i and final time t f , for two constants a and b to set. The minimization problem to compute a and b consists in considering model (2.6)-(2.10) between times t i and t f with known values of S A,(p) (t i ) for all positions passing by the hub between time t i and t f and A ∈ {In, Res, Loa, Out}, and with known functions Φ In,(p) and Φ Out,(p) for every t ∈ [t i , t f ]. We then consider that inner flows are given by (2.17)-(2.19), (2.23)-(2.26) and then a and b must minimize p/T In (p)>t i and T Out (p)<t f

	|ν(p)|,	(4.10)

under the constraint a ≥ 0, b ≥ 0 and dS

Out,(p)

  .1)-(2.27) or (2.30)-(2.45). For this, we consider a network made of n hubs {h 1 , . . . , h n } = H and that the set H is divided into two subsets H mod and H unm (H mod H unm = ∅, H mod H unm = H). Subset H mod contains hubs for which we have on our disposal an instantiation of model (2.1)-(2.27) or (2.30)-(2.45) and H unm contains hubs for which a model of the (2.1)-(2.27) or (2.30)-(

  A i is the set of retailers working with hub #i. The way to connect every instantiation of model (2.1)-(2.27) (indexed by i = 1, . . . , n) consists in translating that goods going out of a given hub towards a second one will later be goods coming from the first hub to the second one. In other words: ) stands for the proportion of mass of goods of position p going from hub #i to #j by travel k l i,j . In particular, Notice that we can consider averaged versions of (4.15) and(4.19)

		Wilfried Despagne & Emmanuel Frénod
	flow of goods going to hub #j. We have:
	Φ	Out,(p) i ∆t (t) = t F j,(p) α∈A i i (s)ds = D α,(p) i L i,j (t) + l=1 ā(k l α∈A i i,j ) t-T l P α,(p) i t-T l i,j +∆t (t) + i,j E	n j=1 i,(p) j	E (s)ds, j,(p) i (t),	(4.14)
		t	∆t	F j i (s)ds =	L i,j l=1	ā(k l i,j )	t-T l i,j +∆t i,j t-T l	E i j (s)ds.
											L i,j
				F	j,(p) i	(t) =	a(k l i,j,p )E	i,(p) j	(t -T l i,j ),	(4.15)
											l=1
	where a(k l i,j,p L i,j
		R λ i (t) =				R	λ,(p) i	(t) , D α i (t) =	D	α,(p) i	(t),
						p∈P os			p∈P os
		B λ i (t) =	p∈P os	B	λ,(p) i	(t) , P α i (t) =	p∈P os	P	α,(p) i	(t),	(4.17)
		F j i (t) =				F	j,(p) i	(t) , E j i (t) =	E	j,(p) i	(t),
						p∈P os			p∈P os
	and,
											n
		Φ In i (t) =				R λ i (t) +	B λ i (t) +	F j i (t),
								λ∈Λ i	λ∈Λ i	j=1 n	(4.18)
		Φ Out i	(t) =				D α i (t) +	P α i (t) +	E j i (t),
								α∈A i	α∈A i	j=1
	and to write,
											L i,j
							F j i (t) =	ā(k l i,j )E i j (t -T l i,j ),	(4.19)
											l=1
	where ā(k l i,j ), satisfying
											L i,j
											j,(p) i	(t) which is the

Concerning the flow of mass of goods of position p leaving "Out Area" of hub #i, we introduce D α,(p) i (t) which is the flow of mass of goods leaving hub #i towards retailer α, P α,(p) i (t) which is the flow of mass of goods directly brought out by retailer α at hub #i, and E where l=1 a(k l i,j,p ) = 1 for any p ∈ P os.

(4.16)

The way to connect every instantiation of model (2.30)-(2.45) indexed by i = 1, . . . , n consists in setting:

l=1 ā(k l i,j ) = 1,

(4.20)

are the proportions of mass of goods going from hub #i to #j via travel k l i,j . ā may be deduced by the history or from data.

Table 12 :

 12 Position parameters over a period of 30 hours.

	Position	Mass	in	Number of	Positions	Unloading	Outbound
	number: p	tons: m(p)	pallets: N p	ahead:	date: T In (p)	date:
					P ub(p)		T Out (p)
	1	5		5	N ull	6	21
	2	8		8	N ull	6	19
	3	2		2	N ull	6	19
	4	3		3	N ull	6	18
	5	10		10	N ull	6	18
	6	5		5	N ull	6.5	20
	7	15		15	N ull	6.5	18
	8	2		2	N ull	6.5	19
	9	3		3	N ull	6.5	18
	10	3		3	N ull	6.5	19
	11	4		4	N ull	7	21
	12	8		8	N ull	7	21
	13	13		13	N ull	7	20
	14	15		15	N ull	7.5	22
	15	3		3	N ull	7.5	22
	16	1		1	N ull	7.5	22
	17	4		4	N ull	7.5	21
	18	14		14	N ull	8	19
	19	7		7	N ull	8	19
	10	6		6	N ull	8	19
	11	2		2	N ull	24	32
	12	6		6	N ull	24	35
	13	10		10	N ull	24	34
	14	5		5	N ull	24.5	33
	15	5		5	N ull	24.5	32
	16	10		10	N ull	25	33
	17	3		3	N ull	25	34
	18	5		5	N ull	25	32
	19	5		5	N ull	26	34
	20	1		1	N ull	26	32
	21	2		2	N ull	26	35
	22	2		2	N ull	26.5	34
	23	5		5	N ull	26.5	34
	24	5		5	N ull	26.5	34
	25	1		1	N ull	27	35
	26	5		5	N ull	27	32
	27	7		7	N ull	27	32
	28	10		10	N ull	27	33
	29	2		2	N ull	27.5	33
	30	5		5	N ull	27.5	32
	31	5		5	N ull	27.5	33
	32	3		3	N ull	28	33
	33	5		5	N ull	28	34
	34	5		5	N ull	28	35
	35	2		2	N ull	28.5


Jean Gregoire

Emilio Frazzoli

Arnaud De

La Fortelle

Tichakorn Wongpiromsarn

This is the supplementary material to the paper: back-pressure traffic signal control with unknown routing rates. It details some model details and proofs not included in the paper due to space limitations. The characterization of the capacity region, the optimality of BP* and the behaviour of the Lyapunov drift under BP control are proved.

I. ROUTING PROCESS ASSUMPTIONS

When a quantity of vehicles arrives at node N a ∈ I(J i ) during slot t, exogenously and endogenously, it is split and added into queues Q ab , b ∈ O(J i ), according to an exogenous routing process R (t) ab , defined for all a, b ∈ N . The arrival process and the routing process are independent, and for all t, R (t)

ab is independent from {Q(τ )} τ ≤t . R (t)
ab takes an integer, returns an integer, and for X ∈ N, b R

ab (X) ≤ X. For all process X(t) such that for all t, R (t) ab is independent from {X(τ )} τ ≤t , there exists a rate r ab ≥ 0 for all a, b ∈ N such that R (t) ab (X(t)) -r ab X(t) is rate convergent with rate 0. As a consequence of the above assumptions: b r ab ≤ 1

II. NETWORK DYNAMICS Consider the network under phase control p(t). Let define the following flow variables:

f ab (t) = min [Q ab (t), µ ab (p(t))]
(2)

f in a (t) = b f ba (t) (3)
The network dynamics under control p(t) is completely described as follows:

Q ab (t + 1) = Q ab (t) + R (t) ab A a (t) + f in a (t) -f ab (t) (4)
Mines ParisTech, Robotics Laboratory, Paris, France † Massachusetts Institute of Technology, Boston, USA * Inria Paris -Rocquencourt, IMARA team, Le Chesnay, France • Ministry of Science and Technology, Thailand

III. CHARACTERIZATION OF THE CAPACITY REGION

A. Stability definition

A key property of queuing systems is stability, defined below: Definition 1 (Stability). The queuing network is stable if each individual queue U satisfies:

lim sup T →+∞ E 1 T T -1 t=0 1 U (t)>V → 0 as V → +∞ (5)
This definition of stability is standard and is applicable to networks with arbitrary inputs and control laws [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF].

B. The capacity region

It is possible to define a capacity region which describes the set of arrivals rates vectors that can be stably handled by the network.

Definition 2 (Capacity region [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF]). Given a routing matrix r, the capacity region Λ r is the closure of the set of all arrival rate vectors λ that can be stabilized by some control.

The following theorem provides a characterization of the capacity region in our particular model: Theorem (Capacity region characterisation). Given a routing matrix r, the capacity region Λ r is the set of arrival vectors λ such that there exists g ∈ Γ satisfying:

∀a, b ∈ N , r ab (λ a + g in a ) ≤ g ab ( 6
)
where Γ is the set of feasible long-term endogenous service rates, defined below:

Γ = Convex_Hull{µ(p) : p ∈ P} (7)
Moreover,

• λ ∈ Λ r is a necessary condition for network stability, considering all possible controls (including those that have perfect knowledge of future events) • λ ∈ int(Λ r ) is a sufficient condition for the network to be stabilized by a control that does not have knowledge of future events.

Proof. The proof is a slightly modified version of the characterization of the capacity region of [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF]. Let Λr denote the set of arrival rates vectors such that there exists g ∈ Γ satisfying:

∀a, b ∈ N , r ab (λ a + g in a ) ≤ g ab (8)
Let prove that λ ∈ Λr is a necessary condition for network stability, considering all possible controls (including those that have perfect knowledge of future events). Suppose that the network is empty at t = 0, using the equation of the dynamics of the network, we obtain:

Q ab (T ) = T -1 t=0 R (t) ab A a (t) + f in a (t) - T -1 t=0 f ab (t) (9)
Suppose the network is stabilized and fix an arbitrary small value > 0. By the network stability necessary condition of [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF], there must exist some finite value V such that at arbitrary large times T , the queues lengths are simultaneously less than V with probability at least 1/2. Hence, there exists a time T such that with probability at least 1/2, the following inequalities hold for all a, b ∈ N :

Q ab (T ) ≤ V (10) V T ≤ (11) T -1 t=0 R (t) ab A a (t) T ≥ r ab λ a - (12)
T -1 t=0 R (t) ab f in a (t) T ≥ r ab f in a (t) - (13) (14)
Define variables g ab =

T -1 t=0 f ab (t)/T . Using the above inequalities together with Equation 9 provides:

r ab λ a + g in a ≤ g ab + 3 (15)
Hence, with probability greater that 1/2, the flows g ab come arbitrary close to satisfying Inequality 8. As a result, there must exist sample paths f ab (t) from which flow variables g ab are defined that satisfy Inequality 15. This proves that λ is a limit point of Λr . As Λr is a compact, it contains its limit points, and we finally obtain: Λ r ⊂ Λr . Now, assume that λ is interior to Λr . Then, as proved in Section IV, one can build a randomized control that stabilizes the network under arrival rates vector λ. Hence, λ ∈ int( Λr ) is a sufficient condition for the network to be stabilized by some control. Combining the two above results proves that Λr = Λ r .

IV. EXISTENCE OF A STABILIZING RANDOMIZED POLICY

The following theorem proves the existence of a stabilizing stationary randomized policy for all arrival rates vectors interior to the capacity region.

Theorem (Existence of a stabilizing randomized policy). Suppose there exists > 0 such that λ + ∈ int(Λ r ), i.e. λ and λ + interior to the capacity region. Then, there exists a stationary randomized control p(t) such that:

E μab (p(t)) -r ab μin a (p(t)) + λ a ≥ r ab a (16)
Proof. For the sake of simplicity, we assume in this proof that E{A a (t)} = λ a and E{R (t) ab (X(t))} = r ab E{X(t)} for all X(t) independent from R (t) ab . It is not true in the general case (it is the particular case of i.i.d. arrivals and routing), and the reader can refer to [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF] for the principle of an extension to the general case using a K-steps Lyapunov drift.

Suppose λ is interior to the capacity region, i.e. there exists a positive vector such that λ + ∈ Λ. By Theorem III-B, there exists g ∈ Γ such that:

∀a, b ∈ N , r ab (λ a + a + g in a ) ≤ g ab (17)
Since g ∈ Γ, it can be expressed as a weighted sum as follows:

g = p∈P w p µ(p) (18)
where weights w p sum to 1. Let define the randomized policy p that selects randomly the phase to apply at every time slot according to probabilities (w p ) p∈P .

It is direct that it will result in a randomized stationary service matrix µ(p(t)) verifying:

E{µ ab (p(t))} = g ab (19)
As a result,

E µ ab (p(t)) -r ab μin a (p(t)) + λ a ≥ r ab a (20)
Now, assume that p(t) is applied to the queuing network. Then, using the equation of the dynamics of the network:

Q ab (t + 1) ≤ max (Q ab (t) -µ ab (p(t)), 0) + R (t)
ab A a (t) + µ in a (p(t)) (21) An inequality holds instead of an equality because the number of vehicles transferred is less or equal to the transmission rate offered by servers.

Squaring both sides and using max 2 (x, 0) ≤ x 2 , we obtain:

Q ab (t + 1) 2 -Q ab (t) 2 ≤ R (t) ab A a (t) + µ in a (p(t)) 2 + µ ab (p(t)) 2 -2Q ab (t) µ ab (p(t)) -R (t) ab A a (t) + µ in a (p(t)) (22)
Define the Lyapunov function

V(Q(t)) = V (t) = a,b Q ab (t) 2 .
Taking expectations, summing over all a, b ∈ N , using independences and noting that E{A} ≤ E{A 2 }, we obtain:

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b Q ab (t)E{µ ab (p(t)) -r ab λ a + µ in a (p(t)) )} (23)
Using Inequality 20, we obtain:

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b r ab a Q ab (t) (24)
Let define η = 2 min a,b r ab a > 0, we finally obtain:

E{V (t + 1) -V (t)|Q(t)} ≤ B -η a,b Q ab (t) (25)
The sufficient condition using Lyapunov drift proved in [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF] enables to conclude stability of the queuing network.

V. OPTIMALITY OF BP* Theorem (Back-pressure optimality). Assuming that pressure functions are linear with strictly positive slopes. Then, BP* is stability-optimal.

Proof. Again, for the sake of simplicity, we assume in this proof that E{A a (t)} = λ a and E{R (t) ab (X(t))} = r ab E{X(t)} for all X(t) independent from R (t) ab . The reader can refer to [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF] for the principle of an extension to the general case using a K-steps Lyapunov drift.

Let θ ab > 0 denote the slope of linear pressure function P ab and Π ab (t) = P ab (Q ab (t)) the evolution of pressures over time. Define the Lyapunov function V(Q(t)) = V (t) = a,b θ ab Q ab (t) 2 and let p(t) denote the control applied to the queuing network. With the same manipulations as for the proof of Theorem III-B, we obtain:

V (t + 1) -V (t) = a,b Π ab (t + 1) 2 -Π ab (t) 2 = a,b θ ab Q ab (t + 1) 2 -Q ab (t) 2 ≤ a,b θ ab R (t) ab A a (t) + µ in a (p(t)) 2 + µ ab (p(t)) 2 -2 a,b θ ab Q ab (t) µ ab (p(t)) -R (t) ab A a (t) + µ in a (p(t)) ≤ B(t)-2 a,b θ ab Q ab (t) µ ab (p(t)) -R (t) ab A a (t) + µ in a (p(t)) (26)
with the upper-bound B(t) defined below:

B(t) = a,b θ ab R (t) ab A a (t) + sup p∈P µ in a (p) 2 + sup p∈P µ ab (p) 2 (27)
Taking expectation and using independences, we get :

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b θ ab Q ab (t)E µ ab (p(t)) -r ab µ in a (p(t))|Q(t) + 2 a,b θ ab Q ab (t)r ab λ a (28)
The upper-bound B is obtained using E{A a } ≤ E{A 2 a }:

B = a,b θ ab r 2 ab A max a + sup p∈P µ in a (p) 2 + θ ab sup p∈P µ ab (p) 2 (29)
By simple sum manipulation, the following identity is obtained:

a,b M ab (g ab -r ab g in a ) = a,b (M ab - c r bc M bc )g ab (30)
Using identity 30, Equation 28 becomes:

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b θ ab Q ab (t) - c r bc θ bc Q bc (t) E {f ab (t)|Q(t)} + 2 a,b θ ab Q ab (t)r ab λ a = B -2 a,b Π ab (t) - c r bc Π bc (t) E {f ab (t)|Q(t)} + 2 a,b Π ab (t)r ab λ a (31)
Now, assume that BP* control p (t) is applied and let V (t) denote the Lyapunov function under p (t). It is assumed that in case of equality when selecting the phase that maximizes the weighted sum, the selected phase p (t)

satisfies µ ab (p (t)) = 0 if W ab (t) = 0.
As a result, we obtain:

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b W ab (t)E {µ ab (p (t)) |Q(t)} + 2 a,b Π ab (t)r ab λ a (32)
By construction of back-pressure control p (t) (see the arg max in the algorithm), p (t) maximizes a,b W ab (t)µ ab (p(t)) over all possible alternative controls p(t). Now, suppose that the arrival rates vector is interior to the capacity region Λ r , i.e. there exists > 0 such that λ + ∈ Λ. Then, as proved in the supplementary material not provided in this paper due to space limitations, there exists a stabilizing stationary randomized control p(t) such that for all a, b ∈ N :

E µ ab (p(t)) -r ab µ in a (p(t)) + λ a ≥ r ab a (33)
Combining the two above statements, taking expectations, and noting that the control p(t) is stationary result in:

a,b W ab (t)E{µ ab (p (t))|Q(t)} ≥ a,b W ab (t)E{µ ab (p(t))} ≥ a,b Π ab (t) - c r bc Π bc (t) E{µ ab (p(t)))} = a,b Π ab (t)E µ ab (p(t)) -r ab µ in a (p(t)) ≥ a,b Π ab (t)(r ab λ a + r ab a ) (34)
Injecting the above result in the Lyapunov drift inequality results in:

E{V (t + 1) -V (t)|Q(t)} ≤ B -2 a,b θ ab r ab a Q ab (t) (35)
Let η = 2 min a,b θ ab r ab a > 0. We finally obtain:

E{V (t + 1) -V (t)|Q(t)} ≤ B -η a,b Q ab (t) (36)
The sufficient condition using Lyapunov drift proved in [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF] enables to conclude stability of the queuing network.

VI. BEHAVIOUR OF THE LYAPUNOV DRIFT IN HEAVY LOAD

CONDITIONS

Theorem (Lyapunov drift under heavy load conditions). Assume λ + ∈ Λ r , BP control is applied and the network is in heavy load conditions, then there exists B, η > 0 such that :

E{V (t + 1) -V (t) | Q(t)} ≤ B -η a Q a (t) (37)
for sufficiently large .

Proof. Again, for the sake of simplicity, we assume in this proof that E{A a (t)} = λ a and E{R (t) ab (X)} = r ab X for all X ∈ N.

Let Π a (t) denote the evolution of P a (Q a (t)) over time and p(t) the control applied to the queuing network. By simple manipulations, we get:

V (t + 1) -V (t) = a θ a (Q a (t + 1) -Q a (t)) 2 + 2 a θ a Q a (t) (Q a (t + 1) -Q a (t)) (38)
As a result, Inequality 38 becomes:

V (t + 1) -V (t) ≤ B(t) -2 a θ a Q a (t) b f ab (t) -R (t) ab A a (t) + f in a (t) (39)
with the upper-bound B(t) defined below:

B(t) = a θ a b R (t) ab A a (t) + sup p∈P µ in a (p) 2 + a θ a b sup p∈P µ ab (p) 2 (40)
Taking expectations and using independences, we obtain :

E {V (t + 1) -V (t) | Q(t)} ≤ B+2 a θ a Q a (t)λ a b r ab -2E    a,b θ a Q a (t) f ab (t) -r ab f in a (t) | Q(t)    (41)
Moreover, by simple sum manipulations, we get the below identity:

E {V (t + 1) -V (t) | Q(t)} ≤ B + 2 a θ a Q a (t)λ a -2 a,b Π a (t) - c r bc Π b (t) E {f ab (t) | Q(t)} (43)
Since c r bc ≤ 1, we obtain:

E {V (t + 1) -V (t) | Q(t)} ≤ B + 2 a θ a Q a (t)λ a -2 a,b (Π a (t) -Π b (t)) E {f ab (t) | Q(t)} (44)
Now, assume that BP control p BP (t) is applied and let V BP (t) denote the Lyapunov function under p BP (t). It is assumed that in case of equality when selecting the phase that maximizes the weighted sum, the selected phase p (t) satisfies µ ab (p (t)) = 0 if W ab (t) = 0. Moreover, by definition of d ab (t), under infinite capacities, f ab (t) = d ab (t)µ ab (p(t)). Hence, we obtain:

E V BP (t + 1) -V BP (t) | Q(t) ≤ B + 2 a θ a Q a (t)λ a -2 a,b W ab (t)E µ ab p BP (t) | Q(t) (45)
By construction of back-pressure control p BP (t) (see arg max in the algorithm), p BP (t) maximizes a,b W ab (t)µ ab (p(t)) over all possible alternative controls p(t). Now, assume λ + ∈ Λ r where is a positive vector. Then, as proved in Section IV, there exists a stabilizing stationary randomized control p(t) such that:

∃g ∈ Γ : ∀a, b ∈ N , E{µ ab (p(t))} = g ab (46)
Combining the two above statements and taking expectations results in:

a,b W ab (t)E{µ ab (p BP (t)) | Q(t)} ≥ a,b W ab (t)E{µ ab (p(t)) | Q(t)} = a,b d ab (t) (Π a (t) -Π b (t)) g ab (47)
Now, assume that the network in heavy load conditions, then d ab (t) = 1 and we obtain:

a,b W ab (t)E{µ ab (p BP (t)) | Q(t)} ≥ a,b (Π a (t) -Π b (t)) g ab (48)
By simple manipulation, we get the following identity: We finally obtain: which can be verified for sufficiently large then η > 0 and the inequality of the sufficient condition for network stability using Lyapunov drift of [START_REF] Neely | Dynamic power allocation and routing for satellite and wireless networks with time varying channels[END_REF] is verified in heavy load conditions. However, it does not imply that the network is stable under BP control since the heavy load assumption is not necessarily verified at any time.

E{V BP (t + 1) -V BP (t) | Q(t)} ≤ B -

P

  a f ab -r ab f in a = a,b P ac r bc P b f ab (42) Using identity 42, Equation 41 becomes:

W-

  -Π b ) g ab = a Π a g out a -g in a ab (t)E{µ ab (p BP (t)) | Q(t)} ≥ a Π a (t) g out a -g in a = a θ a Q a (t) g out a -g in a (50)Moreover, by definition of the input/output flow,g out a -g in a = b g ab -g in a =b g ab -r ab g in a , using the inequalities verified by g ∈ Γ, we obtain:a,b W ab (t)E{µ ab (p BP (t)) | Q(t)} ≥ a θ a Q a (t) b (r ab λ a + r ab a ) -1b r ab g in a (52)Injecting the latter result in Inequality 45 provides:E V BP (t + 1) -V BP (t) | Q(t) ≤ B

Supplementary material to: Back-pressure traffic signal control with unknown routing rates