D Feillet
email: feillet@emse.fr

T Garaix
email: garaix@emse.fr

F Lehuédé

O Péton

D Quadri

A new Consistent Vehicle Routing Problem for the transportation of people with disabilities

Keywords: consistent vehicle routing problem, transportation of people with disabilities, multiple time windows, heuristic branch-and-price

In this paper we address a problem of the transportation of people with disabilities where customers are served on an almost daily basis and expect some consistency in the service. We introduce an original model for the time-consistency of the service, based on so-called time-classes. We then define a new multi-day Vehicle Routing Problem (VRP) that we call the Time-Consistent VRP. We address the solution of this new problem with a Large Neighborhood Search heuristic. Each iteration of the heuristic requires solving a complex vehicle routing problem with multiple time windows and no waiting time which we tackle with a heuristic branch-and-price method. Computational tests are conducted on benchmark sets and modified real-life instances. Results demonstrate the efficiency of the method and highlight the impact of time-consistency on travel costs.

Introduction

Among the diverse categories of problems that have been investigated in operations research, Vehicle Routing Problems probably constitute one of the most emblematic class.

In practice, their purpose is generally to propose optimized distribution plans while delivering goods to customers (be they individuals, firms or retailers, for example). In most of the many models developed for the routing of vehicles, a single period of time is considered (say, a day). This modeling however corresponds to two very different practical situations: either the routes are intended to be performed once and re-computed with new data the next day, or the routes are repeated identically every single day over a long time horizon. However, there is an intermediate situation in which demands are almost identical from day to day. In addition to the usual cost minimization objective, a new criterion might then make sense: limiting the variability of the routes between the different periods of the time horizon. The problem then becomes a multi-day VRP where routes are expected to show some regularity from day to day.

An important illustration of this problem was given by Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF] and concerns fast parcel delivery (involving big companies such as UPS). The ability to assure consistent schedules and to limit the number of drivers that a customer has to deal with can significantly improve the satisfaction of these customers: facilitating receipt operations, improving trust between the delivery companies and their customers, etc. A major advantage can thus be provided in comparison with competing companies.

Another interesting case, which motivated this paper, concerns the transportation of people with mental disabilities from and to day care centers [START_REF] Lehuédé | Un système d'aide à la décision pour planifier les transports vers les établissements médico-sociaux[END_REF]. Because of their lack of autonomy, most of the persons concerned have no vehicle and are not capable of traveling by public transport. Consequently, medico-social centers usually resort to specialized transportation companies to organize and perform the daily trips. Customers need to be transported regularly but sometimes with some small variations during the week. For example, a person may not ask for transportation or may not have the same address every day. In addition, these users are particularly sensitive to changes. The issue of keeping consistent time schedules is therefore fundamental.

In this paper, we investigate how consistency can be considered in vehicle routing problems. Because the expectations of stability are context-dependent, as will be detailed below, we essentially address the transportation of people with disabilities and the stability of the time of services. Since there is a paucity of works on this subject, we first discuss the notion and the modeling of consistency, before introducing a new model (Section 2). We then introduce the Time-Consistent Vehicle Routing Problem and propose a solution approach based on a Large Neighborhood Search (Section 3). Computational experiments evaluating the behavior of this algorithm are then described (Section 4), before the paper is concluded (Section 5).

Modeling the time-consistency of one customer

The issue of constructing consistent vehicle routes is related to many important lines of research in operations research such as robustness or stochasticity [3,[START_REF] Gendreau | A tabu search heuristic for the vehicle routing problem with stochastic demands and customers[END_REF]. However, consistency is generally sought in relation to unexpected events that disrupt the transport.

In contrast, we consider here the consistency of services with respect to deterministic data with slight variations from day to day, which has rarely been investigated in a routing context. Our study mainly compares with the recent paper of Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF] which introduces the consistent vehicle routing problem. Since its first online publication, this problem has been investigated in [START_REF] Kovacs | A template based adaptive large neighborhood search for the consistent vehicle routing problem[END_REF][START_REF] Schneider | Integrated modeling and optimization of the consistent vehicle routing problem with the super node concept[END_REF][START_REF] Tarantilis | Template-based tabu search algorithm for the consistent vehicle routing problem[END_REF]. The consistent vehicle problem clearly shares some similarities with the Courier Delivery Problem (CDP) [START_REF] Sungur | A model and algorithm for the courier delivery problem with uncertainty[END_REF] which, in addition, considers time windows, uncertain service times, uncapacitated vehicles and that frequent customers should not necessarily be served by a single vehicle. A few other studies show some connections with this subject [START_REF] Baniel | Prise en compte d'objectifs de stabilité pour l'optimisation de collecte de déchets[END_REF][START_REF] Benton | The vehicle scheduling problem with intermittent customer demands[END_REF][START_REF] Zhong | Territory planning and vehicle dispatching with driver learning[END_REF]. Finally, the consistent VRP has some relationship with multi-period VRPs [START_REF] Peter | The period vehicle routing problem and its extensions[END_REF], where the assignment of visit days is also considered.

From the customer point of view, consistency improves the quality of service. Having consistent schedules helps the customer to get prepared for the service, to manage material flows, to schedule upstream and downstream tasks, etc. Assigning the same driver to the different visits of a customer facilitates the delivery process and strengthens the relationship between suppliers and customers [START_REF] Day | Improving routing and scheduling decisions at a distributor of industrial gases[END_REF][START_REF] Michael | Assigning delivery routes to drivers under variable customer demands[END_REF][START_REF] Smilowitz | Workforce management in periodic delivery operations[END_REF]. Such criteria can play an essential role in the choice of a supplier [START_REF] Choi | An exploration of supplier selection practices across the supply chain[END_REF].

For drivers, having consistent sequences of visits also presents many advantages. It helps develop a good knowledge of the usual itineraries and traffic conditions, thus enabling drivers to estimate distribution times and potential delays better, reducing transportation times [START_REF] Day | Improving routing and scheduling decisions at a distributor of industrial gases[END_REF][START_REF] Haughton | Measuring and managing the learning requirements of route reoptimization on delivery vehicle drivers[END_REF][START_REF] Smilowitz | Workforce management in periodic delivery operations[END_REF][START_REF] Steinzen | Branching strategies to improve regularity of crew schedules in ex-urban public transit[END_REF]. Driver familiarity with the customer also helps in applying security and administrative procedures [START_REF] Michael | Assigning delivery routes to drivers under variable customer demands[END_REF]. With all these improvements, the transportation process is generally considered more comfortable by the drivers [START_REF] Steinzen | Branching strategies to improve regularity of crew schedules in ex-urban public transit[END_REF].

In the remainder of the paper, we essentially restrict the study to the consistency of schedules. Our assumption is that once a time-consistent solution is computed, finding a reasonably consistent assignment of drivers should be possible. This assumption is checked experimentally in Section 4. The issue of simultaneously optimizing time-and driverconsistency is left for future research.

In order to discuss the modeling of time-consistency, let us consider a customer i, a time horizon D = {1, . . . , |D|} and times of service (or schedule) W = (w 1 , . . . , w |D| ) for customer i over this horizon: the service of customer i starts at time w d on day d. For the sake of simplicity, note that we assume here that customer i is visited every day of the time horizon. Note also that W is assumed to be known. The method to optimize W will be the subject of Section 3.

We first review the model proposed by Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF] and then describe our model which generalizes the former.

Note that we do not consider explicitly customers with several addresses. Actually, time-consistency is important for a given address. A person with two distinct addresses is modeled as two people (one per address).

Modeling time-consistency through time slacks

The time-consistency for customer i is ensured in [START_REF] Groër | The consistent vehicle routing problem[END_REF] by limiting the difference between the latest and earliest service dates. More formally, a time of service variation is measured as:

T C slack (W ) = max 1≤d≤|D| w d -min 1≤d≤|D| w d ( 1
)
This definition is relevant in the context of the small package shipping industry. However, considering other application contexts, the preceding modeling needs to be questioned. A first remark is that it is particularly suitable for small variations in service times. Thus, schedules W 1 = (8:00, 8:03, 8:05, 8:00, 8:02) and W 2 = (8:00, 8:00, 8:05, 8:05, 8:00), such that T C slack (W 1 ) = T C slack (W 2 ) = 5, can be considered equivalent in terms of time-consistency. However, it is less clear for W 3 = (8:00, 9:00, 8:20, 8:30, 8:40) and W 4 = (8:00, 8:02, 8:05, 8:58, 9:00). Although T C slack (W 3 ) = T C slack (W 4 ) = 60, W 4 shows much more regularity, with two sets of almost identical times of service around 8:00 and 9:00. Following this observation, a second important remark is that while W 1 or W 2 undoubtedly define more consistent schedules than W 4 , should the two schedules W 4 and W 5 = (8:00, 8:02, 8:05, 9:58, 10:00), with T C slack (W 4 ) = 60 and T C slack (W 5 ) = 120 respectively, be considered so different? In the following section and in the remainder of this paper, we introduce a new modeling framework for time-consistency which assumes that such schedules are equivalent.

As could be expected, Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF] attain a very high level of consistency in their solutions, thus justifying their modeling. Two factors can be considered determinant in this behavior. First, vehicles deliver 100 to 150 customers everyday. Secondly, many of these customers are visited only once during the time horizon, thus giving many opportunities for substitutions in the routes. As a consequence, consistent schedules can be preserved for the relevant customers, giving rise to a strong time-consistency and mostly avoiding being concerned by the above remarks. From an algorithmic perspective, the chosen approach was to define template routes that respect the so-called precedence principle: "If customers i and j are both served by the same vehicle on a specific day and i is serviced before j, then customer i must receive service before customer j from the same vehicle on all days where they both require service". Note that this principle also deals with driver-consistency.

For the transportation of people with disabilities, the situation is very different. Vehicles only visit a very small number of people, who make transportation requests almost everyday, thus rendering substitution between customers much more complicated and increasing the differences between the service times. This explains why we introduced the modeling developed below.

Modeling consistency through time-classes

From the above discussion, we introduce a new consistency measure that follows these principles: a small gap between service times is not significant; a large gap induces a disturbance, whatever the gap.

Let us consider the parameter δ > 0 reflecting the sensitivity of the customers to schedule variations. In order to evaluate the time-consistency for customer i, we define the notion of time-classes as follows: the times of service w d and w d are allowed to be assigned to the same time-class if |w d -w d | ≤ δ. Hence, the times of service {w 1 , . . . , w |D| } can be divided into time-classes.

It can easily be seen from the above definition that there is no unique way of assigning a set of times of service to time-classes. Let us go back to the example of schedule W 3 = (8:00, 9:00, 8:20, 8:30, 8:40). With parameter δ = 10 one can, for example, construct time-classes {{Mo},{We,Th},{Fr},{Tu}} or {{Mo},{We},{Th,Fr},{Tu}}.

We define the time-consistency T C(W ) of a schedule W = (w 1 , . . . , w |D| ) as the minimal number of time-classes into which w 1 , . . . , w |D| can be ranked:

T C(W ) = min C (2)
subject to

|w d -w d | > δ ⇒ C d = C d (d, d ∈ D), (3)
C d ≤ C (d ∈ D), (4)
C d ∈ {1, . . . , |D|} (d ∈ D), (5)
C ∈ {1, . . . , |D|}. (6)
Constraints (3) impose the assignment of different classes for service times whose difference is greater than δ. Constraints (4) with objective function (2) compute the minimal number of classes required C.

Figure 2 outlines the proposed consistency through time-classes considering schedules W 1 to W 5 introduced in sub-section 2.1, with parameter δ = 10. Note that with this new measure, schedules W 1 and W 2 are still the most consistent, schedules W 4 and W 5 are considered equivalent and schedule W 3 is strongly penalized. The method to compute values T C(W ) and to assign visits to time-classes is presented in the following subsection. Computing the minimal number of time-classes for schedule W then amounts to coloring the vertices in G, such that two vertices connected with an edge are colored with a different color.

Computation of time-classes

While graph coloring problems are generally NP-hard, graph G exhibits a very special structure which implies that the problem is polynomial. This is highlighted by the subsequent equivalent definition of set

E: edge (d, d ) ∈ E if and only if [w d -δ/2, w d + δ/2] ∩ [w d -δ/2, w d + δ/2] = ∅.
As such, the coloring problem can be solved using a greedy algorithm (Algorithm 1). In Algorithm 1, vertices are colored in the increasing order of times of service w d , with the first available color. Let d represent the first vertex colored with the current color. A sufficient condition for a vertex to be colorable using the current color is to be compatible with vertex d.

The time-classes that would be constructed for schedules W 1 to W 5 with δ = 10, applying Algorithm 1, can be observed in Figure 2.

The Time-Consistent VRP

In the remainder of the paper, we first assess how the proposed definition of time-consistency can be handled in a solution approach. Secondly, we evaluate the impact of considering time-consistency on vehicle routes.

The original application includes various specific constraints and criteria that can complicate the readability of the subsequent results. We thus introduce and address a simplified problem. This is called the Time-Consistent Vehicle Routing Problem (TCVRP) and is basically a multi-day VRP with the additional constraint of optimizing time-consistency over days.

An important difference to be noted compared to previous sections is that customers are initially located at the depot and request to be transported back home with regular arrival times of services, with a set of vehicles all departing at the same time (time 0) from the depot. In the real application that motivates this paper, the people with mental disabilities all had to arrive at a given time and expected to be picked up at regular times. From a computational point of view however, the two situations are completely symmetric and equivalent.

We first define the TCVRP in Section 3.1. Then, the proposed solution algorithm is outlined in Section 3.2. Finally, sections 3.3, 3.4 and 3.5 provide details on the different steps of the algorithm.

Definition of the TCVRP

Let us denote D the set of days in the time horizon and consider a complete graph G = (V, A), where V = {0, . . . , n} represents a depot (vertex 0) and the destinations of n customers (vertices 1 to n). Customers need to be transported from the depot to their destination on some days of the time horizon, defined by p d i = 1 if customer i has to be transported on day d, p d i = 0 otherwise. Set K represents a fleet of identical vehicles available at the depot every day of the time horizon. The capacity of a vehicle is limited to Q customers. A travel time t ij > 0 is defined for every arc (i, j) ∈ A. This travel time corresponds equivalently, without loss of generality, to a travel cost. A parameter δ is given to represent the maximal difference of times of service that can be ranked in the same time-class.

Two objectives are pursued. The first is to minimize travel costs, i.e. the sum of the route costs thoughout the horizon. The second is to maximize time-consistency, i.e. to minimize the maximal number of time-classes over the n customers. The number of timeclasses for a customer is defined in equations ( 2)-( 6), where times of service are related to routing decisions.

Note that no time windows are defined for deliveries. In addition, waiting times are not allowed. Without this constraint, waiting could improve time-consistency.

The TCVRP then consists in optimizing these two objectives by finding a set of at most |K| routes every day of the time horizon D such that (i) all the requested deliveries are performed, (ii) all routes satisfy vehicle capacities. Seeing that D is assumed to remain quite small (typically a week), the problem can rather be stated as the problem of solving |D| mono-objective subproblems, where the maximal number of time-classes is limited and successively defined as C = 1, . . . , |D|. The model ( 7)-( 25) presented below shows a mathematical formulation for one of these mono-objective subproblems, with a limited number of time-classes, namely at most C time-classes per customer. Note however that the solution method and experiments presented later deal with the complete (multiobjective) TCVRP. The model is as follows:

Min z = d∈D k∈K (i,j)∈A t ij x d ijk (7) s.t. k∈K j∈V x d ijk = p d i (i ∈ V \ {0}, d ∈ D), (8)
j∈V

x d jik - j∈V x d ijk = 0 (i ∈ V \ {0}, k ∈ K, d ∈ D), (9)
i∈V \{0}

x d 0ik = 1 (k ∈ K, d ∈ D), (10)
i∈V \{0}

x d i0k = 1 (k ∈ K, d ∈ D), (11)
i∈V \{0} j∈V

x d ijk ≤ Q (k ∈ K, d ∈ D) (12)
w d jk ≤ w d ik + t ij + (1 -x d ijk )M ((i, j) ∈ A, k ∈ K, d ∈ D), (13)
w d jk ≥ w d ik + t ij -(1 -x d ijk )M ((i, j) ∈ A, k ∈ K, d ∈ D), (14)
w d i ≤ w d ik + (1 - j∈V x d ijk )M (i ∈ V \ {0}, k ∈ K, d ∈ D), (15)
w d i ≥ w d ik -(1 - j∈V x d ijk )M (i ∈ V \ {0}, k ∈ K, d ∈ D), (16)
w d 0k = 0 (k ∈ K, d ∈ D), (17)
w d 2 i -w d 1 i ≤ δ + M y d 1 d 2 i (i ∈ V \ {0}, d 1 , d 2 ∈ D), (18)
c d 1 i + 1 ≤ c d 2 i + M (1 -y d 1 d 2 i ) (i ∈ V \ {0}, d 1 , d 2 ∈ D), (19)
y d 1 d 2 i ∈ {0, 1} (i ∈ V \ {0}, d 1 , d 2 ∈ D), (20)
c d i ≤ C (i ∈ V, d ∈ D), (21)
x d ijk ∈ {0, 1} ((i, j) ∈ A, k ∈ K, d ∈ D), (22)
w d ik ≥ 0 (i ∈ V \ {0}, k ∈ K, d ∈ D), (23)
w d i ≥ 0 (i ∈ V \ {0}, d ∈ D), (24)
c d i ∈ {1, . . . , C} (i ∈ V, d ∈ D). (25)
The decision variables have the following interpretation: are set to 1 if the times of service of customer i on days d 1 and d 2 are increasing in this order and they are not allowed to belong to the same time-class (w d 2 i -w d 1 i > δ), and 0 otherwise. Constraints [START_REF] Rousseau | A branch-and-price-based large neighborhood search algorithm for the vehicle routing problem with time windows[END_REF] state that customer i is serviced on day d if and only if he or she needs to be transported on that day. Constraints (9) are classical flow conservation constraints. Constraints [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF] and [START_REF] Feillet | New refinements for the solution of vehicle routing problems with branch and price[END_REF] enable vehicles to leave the depot and return to the depot once a day at most, respectively. Constraints [START_REF] Peter | The period vehicle routing problem and its extensions[END_REF] present vehicle capacities. Constraints ( 13) and ( 14) define time variables. They linearize the logical expression

x d ijk : 1 if arc (i,
∀(i, j) ∈ A, ∀k ∈ K, ∀d ∈ D, x d ijk = 1 ⇒ w d jk = w d ik + t ij .
Note that vehicles are not allowed to wait in order to improve time-consistency. Constraints ( 15) and ( 16) define times of service for customers, independently of the vehicle used. They linearize the logical expression

∀i ∈ V \ {0}, ∀k ∈ K, ∀d ∈ D, j∈V x d ijk = 1 ⇒ w d i = w d ik .
Constraints [START_REF] Ioachim | Fleet assignment and routing with schedule synchronization constraints[END_REF] set vehicle starting times at the depot. For every customer, constraints ( 18)-( 20) enforce the assignment of distinct time-classes to service times that differ from more than δ. They linearize the logical expression

∀i ∈ V \ {0}, ∀d 1 , d 2 ∈ D, w d 2 i -w d 1 i > δ ⇒ c d 1 i + 1 ≤ c d 2 i .
Constraints ( 21) limit the maximum number of time-classes allowed for every customer.

Finally, we note that in the case of C ≥ |D|, constraints ( 18)-(20) ( 21) and ( 25) can be removed since every day can be solved independently. Furthermore, the case C = 1 can also be simplified by replacing those constraints by

|w d 2 i -w d 1 i | ≤ δ (i ∈ V \ {0}, d 1 , d 2 ∈ D). ( 26
)
3.2 Solution of the Time-Consistent VRP with a Large Neighborhood Search approach

In this section, we describe the framework of the algorithm proposed for the solution of the TCVRP. Subsequent sections explain the different procedures introduced here.

Time-consistency constraints significantly complicate the problem as they introduce temporal dependences between routes. Very few vehicle routing problems in the literature exhibit similar features. The most closely related problems we were able to find are those with synchronization constraints [START_REF] Bredström | Combined vehicle routing and scheduling with temporal precedence and synchronization constraints[END_REF][START_REF] Ioachim | Fleet assignment and routing with schedule synchronization constraints[END_REF]. As pointed out by Rousseau et al. [START_REF] Rousseau | The synchronized vehicle dispatching problem[END_REF], synchronization constraints considerably limit the applicability of traditional local search techniques. Assessing the feasibility or impact of a very small move can be a large and complex problem. From this observation, we decided to address the solution of the TCVRP with a Large Neighborhood Search (LNS) approach.

The LNS heuristic, first introduced by Shaw [START_REF] Shaw | Using constraint programming and local search methods to solve vehicle routing problems[END_REF] has been successfully applied to vehicle routing problems in recent years [START_REF] Rousseau | A branch-and-price-based large neighborhood search algorithm for the vehicle routing problem with time windows[END_REF][START_REF] Ropke | An adaptive large neighborhood search heuristic for the pickup and delivery problem with time windows[END_REF]. Although our implementation differs from the initial algorithm of Shaw, the proposed algorithm follows the fundamental principle of LNS: the heuristic starts from a feasible solution and alternates between the destruction of a large part of the solution and the reconstruction of the resulting partial solution.

The main principle of the proposed heuristic is to start from a solution that is good from the cost point of view and to reduce the number of time-classes progressively, recording the best found solution for each level of time-consistency (thus solving the bi-objective version of the time-consistent VRP). In order to diversify the search, a perturbation mechanism and a restart procedure are performed periodically.

The initial solution is computed by considering each day independently and solving (heuristically) the corresponding VRPs. The solution obtained tends to minimize the total distance traveled during the time horizon. As time-consistency is completely ignored, it is likely that some customers will have a large number of time-classes.

At every iteration, the destruction phase selects one day for which all routes are destroyed. The reconstruction phase relies on the definition of time constraints that enable the maximal number of time-classes for the customers to be controlled. Our approach is to select a customer whom we attempt to reduce the number of time-classes and to limit the number of time-classes for other customers, to varying degrees depending on their distance from the current maximal number of time-classes. These constraints can be expressed in the form of multiple time windows associated with customers. The problem to solve is then a variant of the VRP that we call the Vehicle Routing Problem with multiple Time Windows and no-wait constraints (VRPmTW-nw). This problem is tackled with a heuristic branch-and-price solution method.

Algorithm 2 outlines the main steps of the LNS heuristic. The LNS algorithm is initialized (Line 1) by solving the VRP defined every day. For this purpose, we apply the heuristic branch-and-price algorithm that is used for reconstruction (see Section 3.4), initialized with routes obtained by applying the Clarke and Wright algorithm. In this case, no time windows are defined for the customers. The main loop (Lines 3 to 13) alternates between the selection of a day to destroy and the reconstruction phase.

Day selection is done at random (Line 5), ensuring however that all days have been selected before one is destroyed again (Line 12). In Line 6, a customer is selected to focus time-class reduction on this customer. The one chosen has a maximal number of timeclasses among customers served on that day. From the customers who satisfy this criterion, we select one for whom the time-class involved on that day has the lowest cardinality. The purpose is to give more chances of deleting that time-class in future iterations.

Following the selection of a day and a customer, an instance of the VRPmTW-nw is defined (Line 7). The way this instance is constructed is presented in Section 3.3. The heuristic branch-and-price method used for the solution of the VRPmTW-nw (Line 8) is detailed in Section 3.4. When possible, the current solution sol is then updated (Line 10) and the LNS iterations are continued until the computing time limit is reached.

When p 1 useless iterations of the LNS operators have not led to any improvement of one of the best known solutions (one per level of time-consistency), diversification or restart components are called for (Line 4). These two components and the selection process between them are presented in Section 3.5. Iterations are counted as useless when they do not succeed in reducing the number of customers with a maximal number of time-classes.

Time-window definition for the VRPmTW-nw

The definition of the instance of the VRPmTW-nw to be solved at every iteration of the algorithm is based on the current solution sol, the day d that was deconstructed and the selected customer i. As explained in Section 3.2, the objective of the deconstruction and reconstruction phases is to reduce the number of time-classes for i and to limit the number of time-classes for other customers according to a set of rules that are described in this section.

In order to describe our scheme, we introduce the following notation. The number of time-classes for j ∈ V \ {0} in sol is denoted T C sol j (D). The number of time-classes for j after the deconstruction of day d is denoted T C sol j (D \ {d}). T C sol (D) is the maximal number of time-classes of the customers during horizon D, that is the time-consistency value for sol. Finally, the construction of the VRPmTW-nw instance relies on a parameter denoted classM ax j . This represents the maximal number of time-classes allowed for j ∈ V \ {0} in the new solution to be constructed.

According to this notation, classM ax j = T C sol j (D) + 1 means that any time of service is possible for customer j on day d (no time constraint). On the contrary, classM ax j = T C sol j (D) generally means that the time of service for customer j on day d has to be compatible with the times of service on other days for this customer. The sole exception is met when T C sol j (D \ {d}) < T C sol j (D), that is when the time of service of customer j on day d in sol was defining a time-class on its own.

In the above paragraph, compatible means that the time of service for customer j on day d can be integrated into one of the existing time-classes for this customer. Given a timeclass whose lowest time of service is t inf and highest time of service is t sup , the condition for the time of service w d j to be compatible with this time-class is t sup -δ ≤ w d j ≤ t inf + δ. This condition can be expressed in the form of a time window condition: w d j ∈ [t sup -δ, t inf + δ]. Being compatible with the different time-classes then comes down to satisfying multiple time windows.

The number and the size of time windows (can be none) assigned to a customer in the travel cost minimization problem on day d depends on how the classM ax parameter has been defined. Let us now describe the proposed definition of classM ax. We recall that i is the customer selected at the current iteration of the LNS algorithm, for whom a decrease in the number of time-classes is expected.

1. Four cases are considered for customers j ∈ V \ {0}, j = i: ii. otherwise classM ax j = T C sol j (D).

(a) if T C sol j (D) = T C sol (D): i. if this condition held in the past p 2 iterations, classM ax j = T C sol j (D) + 1, ii. else, if T C sol j (D) = T C sol (D)

For customer

i, if T C sol (D) > 1, classM ax i = T C sol (D) -1, otherwise classM ax i = 2.
Rules 1(a)i and 1(b)i state that the number of time-classes for customer j can be increased, whereas it remains the same in rules 1(a)ii and 1(b)ii. According to classM ax j , time windows are then defined for j. Two situations are possible:

• classM ax j = T C sol j (D \ {d}) ⇒ [multiple time windows],
• classM ax j > T C sol j (D \ {d}) ⇒ [no time window].

The rationale behind using probability min(p 3 (T C sol (D) -T C sol j (D)), 1) in rules 1(b)i and 1(b)ii is that customers with a large number of time-classes are given less opportunity to increase this number. Parameter p 3 is set to 0 at the beginning of the algorithm and increased by 0.1 at every restart (see Section 3.5) with a maximum value of 0.5. The only possibility that the number of time-classes in the new solution exceeds T C sol (D) comes from rules 1(a)i and 2 (second case), where the purpose is to inject some diversification.

Note finally that the definition of time windows may not actually reduce the number of time-classes of i in the first case of rule 2. In fact, when T C sol i (D \ {d}) = T C sol i (D) = T C sol (D), no solution with T C sol (D) -1 time-classes exists (remember that i is selected from the customers with a maximal number of time-classes on day d). In this case, classM ax i will not exactly represent the maximal number of time-classes for customer i. Definition classM ax i = T C sol (D) -1 will rather be handled by defining one time window per time-class of customer i except for the time-class in which the time of service of customer i on day d was ranked. The purpose is to decrease the size of this time-class and eventually to empty it in further iterations.

Example The preceding principles are illustrated in the following example: we consider a solution for 15 customers on 5 days with δ = 5, which results in a time-consistency T C sol (D) = 3. The selected customer i is customer 7 and day 2 is deconstructed. Columns 2 to 6 of Table 1 present the time schedule of the solution for the 15 customers. Columns T C sol j (D) and T C sol j (D \ {d}) present the number of time-classes of each customer in the current and the deconstructed solutions, respectively. The following columns then present the rule applied to each customer (with an increased number of time-classes allowed for customers 6 and 8) and the time windows defined according to these values ("-" stands for no time window).

Solution of the VRPmTW-nw

The problem to solve at each iteration of the reconstruction phase of the LNS consists in finding a set of routes with minimal cost that serve all customers involved on that day. One additional condition, compared to the standard VRP is that deliveries are restricted by multiple time windows. Also, no waiting time is allowed. We call this problem VRPmTWnw and describe here the solution approach developed for it. While the VRPmTW-nw is new (to the best of our knowledge), a heuristic branch-and-price solution scheme could relatively easily be adapted from existing VRPTW solution schemes.

The principle of the branch-and-price algorithm for the VRPTW [START_REF]GERAD 25th Anniversary Series[END_REF][START_REF] Feillet | A tutorial on column generation and branch-and-price for vehicle routing problems. 4OR[END_REF] is: a set-covering formulation is defined where path variables represent the set of feasible routes. Due to the very large number of possible routes, the formulation is first restricted to a limited subset of variables (restricted master problem) and then extended using a pricing subproblem. At each iteration of the so-called column generation scheme, the linear relaxation of the restricted master problem is solved and the pricing problem is called for, given a vector of optimal dual variables. The aim of the pricing problem is then to determine one or several feasible routes of negative reduced cost. By subtracting the values of dual variables from arc costs, the pricing problem can be expressed as an elementary shortest path problem with resource constraints and solved with dynamic programming. This whole scheme is embedded in a branch-and-price algorithm so as to obtain optimal integer solutions. Branching rules consist in imposing or forbidding arcs in the solution.

Regarding the VRPmTW-nw, most of the previous approaches can be applied in a similar way. The only difference lies in the pricing problem and concerns the feasibility of routes. While extending paths in the dynamic programming algorithm, one has to address multiple time windows and forbidden waiting times. Feasibility is assessed by simply checking that the destination node is reached within one of its time windows. As waiting times are not allowed, the condition that two paths arrive at the destination at exactly the same time is added to the dominance rule. In fact, earlier arrival, which is usually preferred since it gives more opportunity for further extensions, could here prevent some customers from being reached. The consequence is a significantly weakened dominance.

In the process of maintaining elementary paths and reinforcing the dominance rule, label definition includes the computation of the set of reachable nodes [START_REF] Feillet | An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems[END_REF]. Here, reachable nodes are those that were not visited in the partial path represented by the label and whose last time window ends after the earliest possible arrival.

Actually, in view of the way multiple time windows are defined (see Section 3.3), it is possible that a large proportion of customers are not subject to time windows. So as to counterbalance the negative effect of the no-waiting-time condition on the dominance rule, we implement the following mechanism. While computing the set of reachable nodes for a label, we determine that we do not fall into the situation where all the remaining nodes have no time windows. For these labels, the traditional dominance rule is applied. Through these adaptations, the VRPmTW-nw can theoretically be solved exactly. As the branch-and-price algorithm has to be called many times using relatively similar data, a single pool of columns is kept and enriched throughout the successive calls of the algorithm. A filtering algorithm is performed at each new call in order to consider only feasible columns (regarding the current data) in the model. However, the need to solve a VRPmTW-nw instance at every iteration of Algorithm 2 also implies very limited computing times. Hence, a heuristic variant of the previous branch-and-price scheme is developed.

When solving the pricing problem with dynamic programming, the set of labels conserved with every node of the graph is limited to the p 4 = 20 best non-dominated labels, in order of increasing time. Exploration in the dynamic programming algorithm is also limited with a Limited Discrepancy Search policy [START_REF] Feillet | New refinements for the solution of vehicle routing problems with branch and price[END_REF]. Basically, extensions are limited to the p 5 = 2 nearest neighbors (sorted according to the reduced cost). However, extensions to other neighbors are permitted with the limit of at most p 6 = 2 extensions for a complete route.

It is well-known that column generation is subject to slow convergence. In order to avoid long series of column generation iterations without a clear improvement in the objective value, we stop after a given number of iterations with improvements less than p 7 = 1%. This number of iterations is set to p 8 = 5 when a feasible solution has already been found, p 9 = 10 otherwise.

Finally, a simple diving strategy is implemented for the branching phase. When the current optimal solution of the relaxed restricted master problem is fractional, instead of splitting the search into two descendant nodes, we only explore the descendant node where an arc is imposed.

These three heuristic mechanisms avoid, respectively, long pricing problem solution times, a large number of iterations of column generation at each node of the search tree, and a large number of nodes in the search tree, thus ensuring that many instances of the VRPmTW-nw can be solved in a reasonable time. A drawback of these mechanisms is that they involve a relatively large set of parameters (p 4 to p 9 ). The value of the different parameters defined above was set through preliminary experiments, which demonstrated, however, that the method was not very sensitive to these values.

Restart and diversification

The diversification operator works as follows. The first step is to select a customer with a maximal number of time-classes for a maximal number of iterations. This customer is then removed from the solution and reinserted every requested day as the first customer of a route; more precisely, the route that yields the minimal insertion cost among the existing routes on each day. Since vehicles always start from the depot at time 0, the time of service for this customer will always be the same, leading to only one time-class for this customer. This simple operator thus enforces consistency for a customer who was stuck with a large number of time-classes for a long time and starts the search in a completely new direction.

The restart operator replaces diversification when a solution with a time-consistency equal to 1 has been found since the last diversification or restart. It consists in restoring one of the best known solutions as the current solution. Our choice is to restore solutions with any potential number of time-classes. Thus the first call to the restoration operator will generate a solution with |D| time-classes, then |D| -1, etc. When value 1 is reached, the process starts again from |D|. As the LNS algorithm tends to search for solutions with a small time-consistency, this restart operator aims to spend some time improving best known solutions with larger time-consistencies.

Computational experiments

The LNS heuristic was coded in C++ and run on a machine with a 2.5 GHz processor and 3 Go of RAM and evaluated on two types of data instances. First, we considered the so-called small instances of [START_REF] Groër | The consistent vehicle routing problem[END_REF], which comprise 10 instances with 10 or 12 customers over a 3-day horizon. The average level of presence of the customers is around 70%. The instances are named conVRPx-y, with x = 10 or 12 for the number of customers and y = 1, . . . , 5, for the instance index. In these instances, a route duration limit is included, which can easily be handled in the model and the algorithm.

Then, we considered 14 real cases arising from day care centers for mentally disabled children or adults in the area of Nantes (France). The centers have to transport between 15 and 65 people everyday (inbound trips in the morning and outbound trips in the evening). The planning horizon is 5 days. For each original instance, the demand of the transported persons is represented by vectors of 5 binary values: 1 if the person has to be transported on that day and 0 otherwise. We modified the original demand vector so that the level of presence varied from ρ = 0.5 to ρ = 0.9 in steps of 0.1. In order to obtain a level of presence ρ with realistic data, we made the following assumptions: a proportion ρ 2 of the people make a demand everyday; a proportion (1 -ρ 2 ) of the people have a level of presence of ρ 1+ρ ; for the latter category, the level of presence on Tuesdays, Wednesdays and Thursdays is twice that of Mondays and Fridays. The resulting 70 instances are named dataz-x, where z = 10 × ρ and x is the number of customers.

The original instances [START_REF] Lehuédé | Un système d'aide à la décision pour planifier les transports vers les établissements médico-sociaux[END_REF] include additional constraints such as individual time windows or maximum trip duration, several possible addresses during the week, incompatibilities between people, wheelchair constraints, a heterogeneous fleet, etc. We ignored these complicating constraints in order to focus on the consistency aspect.

As explained in the introductory section of the paper, a post-processing phase is carried out to assign drivers to routes in such a way that driver-consistency is optimized. The objective is to minimize the average number of drivers per customer. This is introduced and discussed in [START_REF] Smilowitz | Workforce management in periodic delivery operations[END_REF], where the authors explain that it is representative of customer satisfaction in terms of driver-consistency. lower bounds given by the 2-time-classes problems. All the known optimal solutions were reached at each run of the LNS algorithm, except for instance RconVRP10-5; for this one, our algorithm was not able to find the optimal solution with 2 time-classes; instead solutions of costs 84.5 and 84.69 were found, 8 and 2 times respectively, for an average value of 84.54.

In order to evaluate the results obtained with 1 time-class, Table 3 provides the average and best solutions found by LNS (after the ten runs), the best feasible solution built by the ILP and the best known solution among these solutions ('Best'). Values that match the best known solutions are written in bold. Our heuristic failed to find the best known solution for only two cases. For the RconVRP10-5 instance, it found a sub-optimal solution of cost 96.7 once out of ten runs. For the RconVRP12-2 instance, our heuristic returned a 1-time-class solution with cost 81.34 instead of 81.25 twice out of ten runs.

For three instances out of ten, our heuristic found better results than the ILP formulation. Moreover, since the gap between the best-LNS solutions and the lower boundstaken from the 2-time-classes optimal solutions of Table 2 -is 3.8% on average and 14.5% at maximum, we can be confident in the quality of the solutions computed by our heuristic on the modified small instances.

Comparison with the approach of Groër et al.

Table 4 presents a comparison between our approach and that proposed by Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF]. For this set of experiments the original instances conVRPx-y introduced in [START_REF] Groër | The consistent vehicle routing problem[END_REF] are considered.

The purpose here is not to compare the solution methods but rather to discuss the modeling differences. Indeed, Groër et al. enforce in their model that a single driver is used for every customer, while we do not consider this issue implicitly but rather optimize the number of drivers a posteriori. Also, Groër et al. do not introduce the concept of timeclasses but they consider the maximal difference between the service times. Their heuristic minimizes (implicitly) that difference while it is bounded in their ILP formulation. The chosen limit is 5. This situation corresponds to the set of 1-time-class solutions provided with our heuristic, with parameter δ set to 5.

The optimal costs obtained in [START_REF] Groër | The consistent vehicle routing problem[END_REF] are reported in the column "Opt." of Table 4 1 . The costs of the 1-time-class solutions obtained for these instances with our heuristic are reported in the next column, followed by the observed gap with "Opt.". The difference in cost is not significant on average (3.4%), showing improvements of up to 7.5% when relaxing driver consistency. The impact on driver-consistency can be evaluated by the last three columns of the table, where we report the number of customers with one, two and three drivers, obtained after optimally assigning drivers to the routes in the heuristic solutions. Here, perfect driver-consistency is only reached for instance 'conVRP10-4'. However, the proportion of customers with two drivers never exceeds 30% and no customer has three drivers on these instances. The average losses in driver-consistency (counted in percentage of customers served by more than one driver) equal to 1.8% and 2.6% on instances with 10 and 12 customers, respectively, enable average cost savings of 2.1% and 4.6%. Hence, given the cost reduction obtained, allowing some flexibility in driver-consistency can be worthwhile when possible.

Instance

Opt

Real instances

In this subsection, we investigate the efficiency of our solution approach on the 70 real instances. Remember that in these instances no route duration limit is considered. We first present a set of results (benchmark results) obtained with a computing time of one hour per instance. The objective here is to give some reference results for further research.

Reference results are also provided to evaluate the efficiency of our heuristic with more limited computing times and assess the impact of time-consistency on routing costs.

Benchmark results

Table 5 and 6 details the results obtained by a single run of our heuristic, limited to one hour, on real instances with a time-consistency parameter δ set to 10.

Note that the values given in the first column of this table are the best solutions found with exactly the number of time-classes indicated in the label of the column. In practice, it is clear that a result with fewer time-classes and smaller routing costs would be preferred. For example, on instance data5-15, solutions with 5 and 4 time-classes are actually dominated by the solution with 3 time-classes. Dominated solutions are generally generated by the diversification procedure. When improving cost also decreases the number of time-classes, there may be a large cost difference between the five classes and the four classes solutions. Some dominated solutions may also result from the heuristic solving of the problem or a particular combination of time windows. We decided to present the results as such considering that "non-dominated" solutions can easily be deduced. Statistics on these results and also on "non-dominated" solutions are proposed in Table 7. The first part of the table (Rough results) summarizes the results given in Table 5 and provides the average value for the instances for which a solution was found. In the second part of the table (Non-dominated solutions), values of "dominated" solutions are replaced by values of solutions that "dominate" them (also, the best solution with fewer time-classes is considered when no solution was found). In the last line of the table, gaps indicate the increase in travel costs observed compared to solutions where 5 time-classes are allowed (C = 5).

Several comments can be drawn from Tables 5 to 7. The first rather surprising result is that feasible solutions with a small number of time-classes could be found for all 70 instances. Only some solutions with exactly 5 time-classes were not obtained on the whole set of instances. This latter point is not really surprising since having 5 different classes can sometimes counter-optimize the solution. For example, it is possible that the solution obtained by solving the VRP each day independently (which is optimal in terms of routing) does not introduce any customer with 5 time-classes. In this case, solutions with 5 timeclasses are naturally dominated and not explored during the search.

In addition, we remark that imposing time-consistency does not penalize routing costs so much. As long as 2 time-classes are accepted, the impact on routing costs is less than 1%. When a limit of 1 time-class per customer is strictly enforced, routing costs are only increased by 5.90%, which is not negligible but might be offset by the increase in consistency for the customers.

One can see from Tables 8 and9 that "good" solutions can be obtained relatively quickly using the heuristic (less than ten minutes). After these few minutes, the results are only marginally improved. The only exception concerns solutions with 1 time-class, which are the most difficult to compute and sometimes require many iterations to obtain. The shortest computing times are not assessed as, due to the complexity of the problem to be solved at each iteration of the algorithm, the method is not designed to be efficient in a few seconds.

Impact of parameter δ

Tables 10 and 11 aim to show how solutions are impacted when a tighter gap is fixed. We compare the number of solutions found and the values of the solutions obtained when defining δ = 5 with the previous results (δ = 10). Table 10 is based on rough results, Table 11 considers non-dominated solutions. When computing a gap, only instances for which a solution was found for both values δ are considered. 0.57% 0.66% 0.69% 1.29% 7.09%

Table 11: Gap compared to (non-dominated) solutions found with δ = 10

The main observation that can be drawn from Table 10 is that solutions with one timeclass can still be found most of the time (63 times out of 70). However, we cannot know whether the missing solutions could not be obtained because of our algorithm or they do not exist. Table 11 highlights that when more than 1 time-class is allowed, decreasing δ only slightly impacts the costs. Regarding solutions with 1 time-class, one can observe a clear increase in the costs.

Driver-consistency

Table 12 indicates the average number of drivers assigned to customers. Let us recall that driver-consistency is not explicitly addressed in our model and approach, but is rather handled in post-processing so that the assignment of drivers to the routes computed with the heuristic optimizes the driver-consistency. Table 13 gives the percentage of customers having a given number of drivers out of all the customers of all the solutions (from 5 timeclasses to 1) of the 70 instances. Note that these results are compiled from the solutions presented in Table 5, i.e. the set of rough results. The average number of drivers remains quite limited though not explicitly handled in the algorithm. In addition, we remark that the majority of customers have 1 or 2 drivers.

Conclusion

The content of this paper was initially motivated by the optimization of vehicle routes when transporting people with disabilities to day care centers. So as to capture the specific consistency expectations in this context, a new time-consistency measure, based on time-classes, was defined. From this measure, we introduced the Time-Consistent VRP (TCVRP) and proposed a solution method based on Large Neighborhood Search.

Experiments tend to demonstrate the great difficulty of this new problem and the efficiency of the heuristic proposed. An important result is that, at least for the (realistic) instances used for the computations, very consistent solutions could always be found with a relatively small impact on travel costs. Though driver-consistency was not addressed in the model nor in the solution approach, a fairly good consistency could also be maintained for drivers.

A clear development of this work is to provide efficient exact algorithms to solve the TCVRP. The integer programming formulation proposed in this paper is not solved efficiently by the solver so it would be worth investigating more evolved modeling, valid inequalities and branch-and-cut approaches. From the application perspective, introducing the diverse constraints that were not considered here would also be important. In addition, other possible models of time-consistency could be investigated: for example, it might be interesting to consider the number of time-class changes over successive periods, instead of just counting the number of classes.

Finally, we introduced another new routing problem in the paper, the so-called VRPmTWnw. Given the inherent interest and difficulty of this problem, further research regarding its solution certainly deserves to be conducted.

Figure 1 :

 1 Figure 1: Time-consistency measure of Groër et al.

  It can easily be seen that model (2)-(6) defines a graph coloring problem. Thus, undirected graph G = (D, E) could be introduced, where edge (d, d ) ∈ E if and only if |w d -w d | > δ.

Figure 2 :Algorithm 1

 21 Figure 2: Time-class consistency

Algorithm 2

 2 General scheme of the LNS algorithm 1: sol ← initial feasible solution 2: D ← D 3: repeat 4: every p 1 useless iterations without improvement: diversify or restart 5: select a day d at random in D 6: select a customer i with a maximum number of time-classes among those served on day d 7: define an instance of the VRPmTW-nw for day d according to sol and i 8: s ← solution of the VRPmTW-nw or null if infeasible 9: if s = null then 10: sol ← replace routes of day d in sol by the routes of s if D \ {d} = ∅ then D ← D \ {d} else D ← D 13: until a computing time limit is met

  for no more than p 2 iterations, classM ax j = T C sol j (D), (b) if T C sol j (D) < T C sol (D): i. classM ax j = T C sol j (D)+1 with probability min(p 3 (T C sol (D)-T C sol j (D)), 1),

  drivers 1.43 1.40 1.35 1.30 1.24

Table 1 :

 1 Definition of multiple time windows for a solution with 15 customers, 5 days and δ = 5

	D	1	d=2	3	4	5	T C sol j	T C sol j	Rule	classM ax	Time windows
	V \ {0}						(D) (D \ {d})			for day 2
	1	54.6	54.6	91.8 91.1	91.1	2	2 1(b)ii	2	[44.6;64.6] [81.8;101.1]
	2	12.7				29.4	2	2
	3	23	69			19.1	2	1 1(b)ii	2	-
	4				55.2	55.2	1	1
	5	53.9			7.7		2	2
	6	71.8	15.6	23.8 21.8	21.8	2	2	1(b)i	3	-
	i=7	37.3	94.3	118.3 32.2 117.6	3	2		2 [27.3;42.2] [108.3;127.6]
	8	40.8	40.8	78 77.3	77.3	2	2	1(b)i	3	-
	9	2.4	2.4	2.4 44.4	2.4	2	2 1(b)ii	2	[0;12.4] [34.4;54.4]
	10	84.6		10.6	8.6	8.6	2	2
	11		84.8				1	0 1(b)ii	1	-
	12		29.8	38		51.8	2	2 1(b)ii	2	[28;48] [41.8;61.8]
	13	26.6	26.6	63.8 63.1	63.1	2	2 1(b)ii	2	[16.6;36.6] [53.8;73.1]
	14		108.6		17.9		2	1 1(b)ii	2	-
	15			8.1			1	1

Table 3 :

 3 They also propose several integer programming 1-time-class solutions on modified small instances (δ = 5)

		Avg LNS Best LNS	ILP	Best
	RconVRP10-1	92.91	92.91	92.91	92.91
	RconVRP10-2	80.99	80.99	82.83	80.99
	RconVRP10-3	94.37	94.37	94.37	94.37
	RconVRP10-4	94.09	94.09	94.53	94.09
	RconVRP10-5	96.08	96.01	96.01	96.01
	RconVRP12-1 104.40	104.40	106.87	104.40
	RconVRP12-2	81.27	81.25	81.25	81.25
	RconVRP12-3	83.12	83.12	83.12	83.12
	RconVRP12-4 101.91	101.91	101.91 101.91
	RconVRP12-5	89.25	89.25	89.25	89.25

Table 4 :

 4 Heuristic solutions and driver-consistency on small instances (δ = 5)

	.(*) LNS		Number of customers with
	Cost	Cost	Gap 1 driver 2 drivers 3 drivers
	conVRP10-1 122.03 118.70 2.8%	7	3	0
	conVRP10-2 99.07	96.57 2.6%	8	2	0
	conVRP10-3 123.41 123.16 0.2%	9	1	0
	conVRP10-4 126.89 126.89 0.0%	10	0	0
	conVRP10-5 109.31 104.04 5.1%	7	3	0
	conVRP12-1 144.02 143.23 0.6%	10	2	0
	conVRP12-2 89.54	83.27 7.5%	9	3	0
	conVRP12-3 119.69 113.74 5.2%	9	3	0
	conVRP12-4 141.37 137.26 3.0%	10	2	0
	conVRP12-5 117.42 110.17 6.6%	9	3	0

Table 12 :

 12 Average number of drivers per customer

	Number of drivers	1	2	3	4	5
	Percentage of customers 69.06% 25.9% 4.56% 0.46% 0.00%

Table 13 :

 13 Percentage of customers having a given number of drivers

Optimal values given in[START_REF] Groër | The consistent vehicle routing problem[END_REF] include service durations. They are discarded in the costs of Table4.

Acknowledgements

This research was funded by the GdR RO (Operations Research Group of the Scientific Research National Center).

formulations. In our computations, we used their (DC) model which minimizes the average number of drivers transporting each customer.

Small instances

The aim of this subsection is threefold: (i) to show that the problem is difficult to solve in practice (i.e., the limits of the integer linear model); (ii) to validate the proposed heuristic in comparison with the best solutions found; (iii) to compare the proposed model with Groër et al. [START_REF] Groër | The consistent vehicle routing problem[END_REF].

Comparison with the integer linear program

Table 2 first evaluates the ability of the integer linear programming formulation introduced in Section 3.1 to solve the problem. As it seems almost incapable of solving any of the instances, we temporarily relax vehicle capacity constraints, maximum route length constraints and assume unitary demands. These instances, derived from the original ones, are named "RconVRPx-y". The sensitivity parameter δ is set to 5. Three runs of 10 hours using ILOG-CPLEX 12.0 are performed on each instance, setting the maximal number of time-classes C -in inequality ( 21) -to 3, 2 and 1, respectively. The optimal solutions found are presented in the "ILP" section of the table. Optimal solutions are compared to the solutions provided by our heuristic in 5 minutes of computing time (each run provides a solution for every number of time-classes): the "Avg LNS" part presents the average results over ten runs of the heuristic.

ILP

Avg LNS Time-classes: Regarding the solution method, it is interesting to add that an average of 6000 iterations could be performed in one hour, thus allowing a large exploration of the solution space (remember that in one iteration, one fifth (one day) of the solution can be completely redrawn).

Efficiency of the algorithm with shorter computing times

In Tables 8 and9, we compare the number of solutions found and the values of the solutions, for computing times varying between 60 seconds, 600 seconds and one hour respectively.


A Muller

C Pontonnier

G Dumont

MusIC method enhancement by a sensitivity study of its performance: application to a lower limbs musculoskeletal model

Keywords: Musculoskeletal analysis, Motion, Computation time, Database, Biomechanical feedback, Interpolation

published or not. The documents may come

MusIC method enhancement by a sensitivity study of its performance: application to a lower limbs musculoskeletal model

Antoine Muller, Charles Pontonnier, Georges

Introduction

Providing a feedback based on biomechanical quantities during experimental sessions is a real outcome for rehabilitation, ergonomics or training applications. The interaction between the subject, the experimenter (e.g. trainer, ergonomist or clinician) and the biomechanical data (e.g. muscle forces) in experimental sessions is a promising way to adapt and enhance the subject performance according to the results of the analysis. However, such applications imply a fast computation of the biomechanical quantities to be observed -particularly muscle forces. The main challenge in estimating muscle forces is to treat the muscular redundancy question that leads to an underdetermined problem. This problem is classically solved by gradient-based numerical optimization schemes. In such an optimization scheme, the minimized cost function aims at mimicking the central nervous system behavior [START_REF] An | Determination of muscle and joint forces: a new technique to solve the indeterminate problem[END_REF][START_REF] Crowninshield | A physiologically based criterion of muscle force prediction in locomotion[END_REF][START_REF] Rasmussen | Muscle recruitment by the min/max criterion -A comparative numerical study[END_REF]) and the constraints represent the dynamical equilibrium of the musculoskeletal model and the maximal forces that the muscles can develop [START_REF] Challis | An analytical examination of muscle force estimations using optimization techniques[END_REF][START_REF] Challis | Producing physiologically realistic individual muscle force estimations by imposing constraints when using optimization techniques[END_REF][START_REF] Erdemir | Model-based estimation of muscle forces exerted during movements[END_REF]. Numerical optimization under non-linear constraints is usually computationally expensive. Thus, in most cases, the analysis is computed off-line thanks to pre-recorded and pre-processed motion data, that is unadapted to provide a feedback during the experimental session.

Several studies proposed solutions to provide a real-time visualization of muscle forces generating a given gesture. [START_REF] Murai | Musculoskeletal-see-through mirror: Computational modeling and algorithm for whole-body muscle activity visualization in real time[END_REF] proposed an adapted optimization scheme using electromyography (EMG) data and gathering muscles defined in the musculoskeletal model into functional groups in order to reduce the number of unknowns to optimize. This simulation required the use of EMG-sensors and made an important assumption by computing a unique activation level for all the muscles of a given functional group. As an alternative, ( [START_REF] Van Den Bogert | A realtime system for biomechanical analysis of human movement and muscle function[END_REF] proposed to use a neural network trained to mimic quadratic optimization. The computation time was limited to ensure real-time computation. Both of these works proposed real-time musculoskeletal forces estimation thanks to the methods and assumptions described above, that limited the accuracy of their results.

To limit the computation time, (Pontonnier andDumont 2009, 2010) proposed to interpolate pre-computed solutions from a muscle forces database. Limitations remained numerous since the method was applied to a single joint with only 4 muscles and a large work had to be done to extend the method to musculoskeletal models with more joints, muscles and muscular coupling -i.e. multi-articular muscles. Based on the same idea, [START_REF] Muller | The MusIC method: a fast and quasi-optimal solution to the muscle forces estimation problem[END_REF] proposed the MusIC (Muscle forces Interpolation and Correction) method as a global solution to the muscle forces estimation problem. The results were very encouraging since the on-line computation frequency (number of samples computed per second) was multiplied by 10 with a mean error of 4% with regard to the optimization scheme. Three different cost functions (e.g. min-max, quadratic, cubic) were tested without any restriction.

Several issues needed to be solved to make the MusIC method being a credible solution to be used in a feedback setup. First, the method was applied on a 2 degrees of freedom (DOF) planar representation of the human arm (shoulder and elbow flexion) actuated by 12 muscles. This model was useful to evaluate the MusIC method performance for a large number of configurations but was not representative of a real musculoskeletal analysis, asking for more accurate and 3-dimensional models. Second, the off-line computation time -corresponding to the database generation time -was about 1 hour. Since one database has to be generated per subject, as each scaled model needs its own muscle force database in this method, the database generation has to be computed after the calibration of the model, that is in every case the first step of the experimental session. Therefore, 1 hour of waiting time per subject is not realistic for an effective use of the method in an experimental session with direct biomechanical feedback. Moreover, the database generation time grows with the complexity of the models to simulate due to the number of lines of action and degrees of freedom.

In order to circumvent these limitations, the current paper proposes to adapt and enhance the usability of the MusIC method for a complex 3-dimensional musculoskeletal model. First, the MusIC method has been applied to a lower limbs musculoskeletal model available in the AnyBody Modelling System [START_REF] Damsgaard | Analysis of musculoskeletal systems in the AnyBody Modeling System[END_REF] to assess both its accuracy and computation time. Next, a sensitivity study has been conducted to find a trade-off between the MusIC method accuracy and the off-line computation time needed for the database generation. This method is able to find a priori -without generating databases -the databases exhibiting the best trade-offs between accuracy and computation time.

The article begins with a short presentation of the MusIC method. Second, assumptions are formulated to test if a database density can be decided without generating the associated database. Third, additional parameters are introduced to adjust the database structure in order to better control the database discretization. Then, the method is applied on a 3-dimensional lower limbs musculoskeletal model on a large set of database densities and evaluated by generating results of the MusIC method for all of these densities. Last, the results of accuracy and computation time are presented and discussed.

Materials and methods

In the following section, we consider a generic musculoskeletal model composed of n b bodies, articulated by n j joints and actuated by n m muscles. The (n j × 1) vector q contains the joint coordinates, the (n j × 1) vector Γ contains the joint torques and the (n m × 1) vector F contains the muscle forces.

The MusIC method

From joint coordinates q and joint torques Γ, the MusIC method aims at computing a muscle forces vector solution F on the basis of a precomputed database. The method is composed of two steps: the off-line database generation -only based on the musculoskeletal model -and the on-line muscle forces estimation -based on the database and on motion data. The database generation consists in storing solutions of the optimization problem for a set of joint configurations. From a current configuration, a first estimation is computed thanks to a linear interpolation in the database. Then, a constrained optimization corrects it to satisfy the dynamic equilibrium -respect of equations linking the joint torques and the muscle forces -and the physiological limits of the muscles -muscles can only pull with maximal forces. Readers are referred to [START_REF] Muller | The MusIC method: a fast and quasi-optimal solution to the muscle forces estimation problem[END_REF] for more details about the MusIC method.

Influence of database density

The influence of the density of the database on the accuracy and the off-line computation time of the MusIC method is investigated in this part. For a given joint configuration, the MusIC method uses the closest configurations stored in the database in order to interpolate the related muscle forces. The more the joint configuration is discretized in the database, the more the result given by the interpolation is close to the optimization one. Moment arms are variables which influence the results of the optimization. So we can assume that investigating the error made on moment arms computation throughout the joint space discretization may explain most of the error provided by the interpolation method. To assess the database accuracy and computation time without generating them, we assumed the two following hypotheses to be discussed:

[H1] the accuracy of the MusIC method is directly correlated to the accuracy of the moment arms computation throughout the joint space discretization;

[H2] the densities with highest trade-off between accuracy and off-line computation time for a scaled model are the same as whose defined on a generic model. We define the term generic model as a non-scaled model. The parameters of the scaled model are calibrated according to the studied subject.

To investigate these assumptions, we propose to use a generic musculoskeletal model and to evaluate the error made on moment arms computation during the joint space discretization. Therefore if hypotheses [H1] and [H2] are valid, this error should be correlated to the error made in the MusIC method with a scaled model. Thus, we propose to validate the assumptions by evaluating the accuracy of the muscle forces estimation with regard to the discretization of the database for a given set of motions and subjects.

Database structure adjustments

Let us consider a joint to be analyzed, called primary joint. The secondary joints of this primary joint are those which have an impact on the primary joint actuation, i.e., in our case, which have an influence on the moment arms values of the muscles actuating the primary joint. Figure 1 illustrates the primary and secondary joints concepts on a simple musculoskeletal model.

C(1) = 2 C(2) = 1, 3 C(3) = 2 C(4) = ∅ Figure 1.
Secondary joints associated to each primary joint of a simple musculoskeletal model. For each primary joint, the secondary ones are those which have an influence on moment arms values of muscles actuating the primary joint. C represents the index function of secondary joints. On this example, a bi-articular muscle between joint 1 and joint 2 induces that joint 2 is the secondary joint of joint 1. The bi-articular muscles between joint 1 and joint 2 and between joint 2 and joint 3 induces that joints 1 and 3 are secondary joints of joint 2. The bi-articular muscle between joint 2 and 3 induces that joint 2 is the secondary joint of joint 3. Since all the muscles actuating joint 4 are mono-articular, no secondary joint is associated to this joint.

The precomputed database of the MusIC method contains a set of solutions of the optimization problem. These solutions are solved by considering each joint separatelyas primary joints, therefore one sub-database is generated per joint (left part Figure 2). In a sub-database, only the primary joint and its secondary ones are considered. The accessible, constrained by the joint limits, joint space associated to these joints is discretized to evaluate the evolution of the muscle forces estimation solution according to the joint configuration (central part Figure 2). For each of these joint configurations, the solution of the optimization problem is stored (right part Figure 2).

Database

Sub-database i

Discretization of joint space including the primary joint i and the secondary joints

Configuration j

Solution of the optimization problem

Sub-database n j

Sub-database i

Sub-database 2 Sub-database 1

One sub-database associated with each joint

q i-1 q i+1
Secondary joints (k 2 intervals)

q i
Primary joint (k 1 intervals) We estimate the database density thanks to the number of samples that it contains. For each sub-database, we assume that the angular position of the primary joint has more impact on the moment arms values of this joint than the angular position of the secondary joints. As such, we propose to use different sampling density to discretize the primary joint amplitude and the secondary ones. The database density is defined thanks to the parameters (k 1 , k 2 ), where k 1 is the number of intervals for the discretization of the primary joint space and k 2 is the number of intervals for the discretization of the secondary ones. A regular sampling will be used in the following.

Databases accuracy evaluation

We evaluate the error made on moment arms computation in each database by interpolating moment arms for a large set of joint test configurations. Each joint coordinate is discretized into k t testing values between its physiological limits that leads to k nj t joint test configurations. For each joint database (defined by k 1 and k 2 values), moment arms are linearly interpolated from the joint configurations of the database for the k nj t joint test configurations. For the j-th configuration, the interpolated value of the moment arms R interp,j is computed and compared to its reference value R ref,j . R ref,j is computed thanks to its analytical expression, obtained from the analytical description of muscles anatomical references (origins, insertions and via points).

Finally, for a specified (k 1 , k 2 ) discretization, the error made on moment arms computation is evaluated by averaging the absolute difference between R ref and R interp for all k nj t joint configurations (1). Thus, with regard to [H1], the lower this difference is, the higher the accuracy of the MusIC method should be.

∆R dis = 1 k nj t k n j t j=1 |R ref,j -R interp,j | (1)
∆R dis is defined as the relative error of ∆R dis with regard to the minimal value of all databases discretizations (2). ∆R min dis corresponds to the minimal value of ∆R dis . ∆R dis is called the moment arm error.

∆R dis = ∆R dis ∆R min dis -1 (2)
Once the moment arms errors are evaluated, [H1] and [H2] are challenged through the application of experimental data to several databases generated with different densities. For each database, the accuracy of the MusIC method is associated with the off-line computation time (time used to generate the database).

As the MusIC method aims at emulating an optimization method, we evaluate its accuracy by computing, with the interpolated muscle forces solution, the value of the cost function cf of the optimization method used to generate the solutions of the database. To facilitate the results exploitation, as done previously for the moment arms computation, ˜ cf is defined as the relative error on cf with regard to the minimal cf value of all the databases of the study (3). min cf corresponds to the minimal value of cf . ˜ cf is called the cost function error.

˜ cf = cf min cf -1 (3)
To validate the assumptions [H1] and [H2], a Spearman correlation coefficient between moment arm errors ∆R dis and cost function errors ˜ cf for a given database discretization is computed.

Last, the database generation time is stored for each database generation. Indeed, the accuracy of the MusIC method is assumed to depend on the database density ([H1]). In a similar manner, since each sample computation in the database needs to solve an optimization problem, the database generation time -off-line computation time -grows with the number of samples to populate the database and would normally be correlated with the accuracy of the method.

Application

The method proposed above has been applied on a musculoskeletal model commonly used for analysis of lower limbs activities (Figure 3) -e.g. in the AnyBody Modelling System software [START_REF] Damsgaard | Analysis of musculoskeletal systems in the AnyBody Modeling System[END_REF]. It is composed of 82 muscles and actuates 6 joints containing 12 DOF -3 for the hip, 1 for the knee and 2 for the ankle. This model uses a majority of poly-articular muscles and allows the method to be tested with a significant muscular coupling.

First, a generic version of this model close the 50th percentile -a size of 180 cm was used to assess the error realized on the interpolated moment arms. From this generic model, a set of 28 joint space discretizations was used. Each discretization corresponds to a couple of (k 1 , k 2 ) values -(2, 2), (3, 2), (4, 2), (5, 2), (8, 2), (10, 2), (15, 2), (20, 2), (30, 2), (3, 3), (4, 3), (5, 3), (8, 3), (10, 3), (15, 3), (20, 3), (30, 3), (4, 4), (5, 4), (8, 4), (10, 4), (15, 4), (20, 4), (5, 5), (8, 5), (10, 5), (15, 5) and (20, 5). For each discretization, the density was stored and the moment arm error was evaluated thanks to the equations ( 1) and (2). For this application, the parameter k t was set at the maximal value of k 1 and k 2 : k t = 30. As the generic musculoskeletal model used is symmetric, the moment arm error is computed by only considering the muscles of one leg.

Number Muscles 82

Mono-articular 8

Bi-articular 16

Tri-articular 40

Quadri-articular 18 Second, the hypotheses were challenged by applying the MusIC method to scaled models using experimental data. To this end, 10 male participants (27.5±5 years old, height: 180±7 cm, mass: 74.1±7.3 kg) participated to an experimentation. All subjects signed an informed consent form before participation and an anonymization protocol was followed for data safeguarding. Each subject was asked to mimic an avatar performing a standardized motion in a video (Figure 4). This video guide was designed to activate sequentially each degree of freedom. The motion was captured thanks to 16 Vicon R cameras (100 Hz) and 47 motion capture markers placed on standardized anatomical landmarks as recommended by the International Society of Biomechanics (ISB) [START_REF] Wu | Isb recommendation on definitions of joint coordinate system of various joints for the reporting of human joint motion-part I: ankle, hip, and spine[END_REF][START_REF] Wu | Isb recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion-part II: shoulder, elbow, wrist and hand[END_REF]. Moreover, the ground reaction forces were measured with two AMTI R force platforms (1000 Hz), each supporting one foot.

From experimental markers position, the geometrical parameters of the model were scaled thanks to an optimization-based calibration method [START_REF] Muller | A simple method to calibrate kinematical invariants: application to overhead throwing[END_REF]. The inertial parameters of the model were deducted and scaled from an anthropometric rule [START_REF] Dumas | Adjustments to McConville et al. and Young et al. body segment inertial parameters[END_REF]). Each scaled model corresponds to this calibrated model. The joint coordinates were computed using a global optimization method [START_REF] Lu | Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints[END_REF] and the joint torques were computed using a recursive Newton-Euler algorithm [START_REF] Featherstone | Rigid body dynamics algorithms[END_REF] where measured ground reaction forces were considered as external forces.

For each subject, the 28 datab ses, corresponding to the (k 1 , k 2 ) couples described above, were generated only for validation. The corresponding time generation was stored for each database. The cost function used to generate the databases was a classical polynomial criterion corresponding to the cubic sum of the normalized muscle forces. Muscles forces were normalized by their maximum values based on anthropometric tables (no force-length nor force-velocity relationships were used here). Thanks to the MusIC method [START_REF] Muller | The MusIC method: a fast and quasi-optimal solution to the muscle forces estimation problem[END_REF], the muscle forces were estimated from each database previously described. The accuracy was then evaluated thanks to equation (3).

Results and discussion

Optimal densities

Considering the generic model, the moment arm error of a given database -evaluated by ∆R dis -has been plotted with regard to the number of samples of the associated database, directly computed from k 1 and k 2 (Figure 5). As expected, the moment arm error is directly linked to the database density. The higher the level of discretization is, the higher is the accuracy of the moment arms interpolation. Densities with k 2 = 2 belonging the Pareto frontier contain a small number of samples. Therefore, the corresponding off-line computation time is low. However, the errors made on moment arms computation can be more than ten times higher than the ones obtained the most accurate densities. In opposite, optimal densities with k 2 = 5 enable to obtain the highest accuracies. In return, with regard to the number of samples, the off-line computation time is higher. The optimal densities with k 2 = 3 or k 2 = 4 provide an interesting trade-off between the accuracy of the moment arms evaluation and the number of samples. Thus, a Pareto frontier appears in Figure 5.

Figure 6 depicts the MusIC method results obtained from the experimental data. First of all, the generation of the databases shows that the time generation is directly linked to the number of samples.

Then, considering figures 5 and 6, the moment arm errors ∆R dis and cost function errors ˜ cf exhibit similar trends with regard to the database discretization and computation time. The Spearman correlation between the mean of ˜ cf and ∆R dis is 89% for each database. These results seem to confirm assumptions [H1] and [H2] previously described. Indeed, these strong correlations indicate that the MusIC method accuracy is correlated to the moment arms interpolation accuracy. Moreover, the Mu-sIC method accuracy is barely affected by the model scaling since this correlation is computed between moment arms interpolated from a generic model and the MusIC method results applied to scaled models. In addition, the Spearman correlation between the standard deviation of ˜ cf and ∆R dis is 79% for each database. Thus, more 8 10 5 0 0 15,000 30,000 45,000

Number of samples of the discretization

∆R dis k 2 = 2 k 2 = 3 k 2 = 4 k 2 = 5
Pareto frontier we have chosen to associate one color for all (k 1 , k 2 ) parameters with the same k 2 value. Thus, the blues crosses correspond to the (2, 2), (3, 2), (4, 2), (5, 2), (8, 2), (10, 2), (15, 2), (20, 2) and (30, 2) values; the red stars correspond to the (3, 3), (4, 3), (5, 3), (8, 3), (10, 3), (15, 3), (20, 3) and (30, 3) values; the green dots correspond to (4, 4), (5, 4), (8, 4), (10, 4), (15, 4) and (20, 4) values and the black squares correspond to (5, 5), (8, 5), (10, 5), (15, 5) and (20, 5) values. The number of samples of the discretization is directly related to k 1 and k 2 .

the density is important, more the magnitude of the error is low. The validation of the two assumptions [H1] and [H2] indicates that optimal densities of the database can be directly deducted from the moment arm error. Since assumptions [H1] and [H2] are validated, the figure 5 can be used to find an optimal database density (k 1 , k 2 ) a priori, meaning there is no need to compute the database to have an realistic idea of its accuracy. Therefore, depending on the target application, one may choose a given density for a given level of performance (computation time, accuracy). For example, for the model used in the current study, a cost function error of 0.1% could be obtained with a database generation time lower than 10 minutes by using the database with the (4, 4) density.

Practical considerations

The application presented in the current paper has been computed on a classical laptop (Intel Core i7 2.8 GHz processor) with Matlab R . For this musculoskeletal model (exhibiting 82 muscles, 6 joints and 12 DOF), once the database generation is done, the on-line muscle forces computation frequency is about 58Hz with the MusIC method and 1.6Hz with a classical optimization method -for comparison purpose, the stopping criterion of the optimization method is selected to obtain the same level of accuracy as the MusIC method using the highest database density. Therefore, with a sufficiently low database generation time, the MusIC method may be considered as a feedback tool for users or experts in training or ergonomics applications since it can provide in a really fast time a reliable information about the muscle forces involved in the studied motion.

The use of the cost function error as a metric to test [H1] does not guarantee the same muscle forces repartition. However, as shown in [START_REF] Muller | The MusIC method: a fast and quasi-optimal solution to the muscle forces estimation problem[END_REF]), a good Database generation time (min)

Figure 6. Application of the MusIC method to the experimental data. Cost function error ˜ cf for each database discretization are plotted according to the corresponding database generation time. Let us notice that the vertical axis is plotted on a logarithmic scale. The blues box plots correspond to the (2, 2), (3, 2), (4, 2), (5, 2), (8, 2), (10, 2), (15, 2), (20, 2) and (30, 2) values; the red ones correspond to the (3, 3), (4, 3), (5, 3), (8, 3), (10, 3), (15, 3), (20, 3) and (30, 3) values; the green ones correspond to (4, 4), (5, 4), (8, 4), (10, 4), (15,4) and (20,4) values and the black ones correspond to (5, 5), (8, 5), (10, 5), (15, 5) and (20, 5) values.

optimality of the cost function tends to obtain muscle forces close to those computed thanks to an optimization method.

When performing a musculoskeletal analysis session based on motion capture, the MusIC method can be used only if a database has been generated on the basis of the scaled model. Thus, this database generation has to be done after the scaling of the model. Since this scaling implies an experimental data collection [START_REF] Muller | Non-invasive techniques for musculoskeletal model calibration[END_REF], the decrease of the generation time allows the database to be generated on the fly during an experimental session, between the experimental calibration of the model and the capture of the motions to be analyzed. This computation time might be used by the experimenters to explain the experimentation protocol to the subject. In this context, the database with the density corresponding to (k 1 , k 2 ) = (4, 4) seems to be an good trade-off. It allows to obtain a cost function error of about 0.1% with a database generation time lower than 10 minutes.

The influence of the database density on the on-line computation time -evaluated thanks to the computation frequency -has not been mentioned yet. In the current application, the average computation frequency for all the motions, subjects and databases is 58Hz(±3Hz), tested on the same laptop. Thus, with regard to this small standard deviation, we can consider that it does not influence the choices made on database densities.

Limits and perspectives

Optimal densities obtained in this study are specific to this musculoskeletal model. However, the method proposed in this paper allows to find the optimal densities without generating any database. For this, no experimental data is necessary. Thus, this study could be performed, a priori to any experiment, and included in the experimental design and in the model choice.

Moreover, for a given joint space discretization, the accuracy of the MusIC method applied to a scaled model is correlated to the accuracy of the moment arms interpolation applied to a generic model. The sizes of the subjects range between 171 cm and 190 cm. The differences between the scaled and the generic model are restricted, since the generic model was representing a subject of 180 cm. The validity of the assumption [H2] has to be tested with other more distant further from the 50th percentile than these tested subjects. Moreover, the assumption is only confirmed for a geometrical optimization-based scaling and for an anthropometric inertial scaling. In order to enhance the results of the musculoskeletal study, other works suggest to scale the inertial parameters from motion and ground reaction forces [START_REF] Ayusawa | Real-time implementation of physically consistent identification of human body segments[END_REF][START_REF] Hansen | An individual and dynamic body segment inertial parameter validation method using ground reaction forces[END_REF][START_REF] Jovic | Humanoid and human inertia parameter identification using hierarchical optimization[END_REF] and the muscular parameters from isokinetic ergometer measures [START_REF] Garner | Estimation of musculotendon properties in the human upper limb[END_REF][START_REF] Van Campen | A new method for estimating subject-specific muscle-tendon parameters of the knee joint actuators: a simulation study[END_REF][START_REF] Heinen | Prediction of muscletendon parameters based on isokinetic measurements[END_REF][START_REF] Heinen | Muscle-tendon unit scaling methods of hill-type musculoskeletal models: An overview[END_REF]. The use of a complete scaling -geometrical, inertial and muscular parameters [START_REF] Muller | Non-invasive techniques for musculoskeletal model calibration[END_REF]) -is therefore an interesting complementary step to validate assumption [H2], and also to increase the accuracy of the musculoskeletal analysis.

The muscle model used in this paper is simplified since the maximal available forces are considered as constants. Therefore the muscles forces estimation was only influenced by the joint moment arms. A force-length or force-velocity relationships could be used to model the muscle physiological behavior [START_REF] Zajac | Muscle and tendon: properties, models, scaling, and application to biomechanics and motor control[END_REF][START_REF] Rengifo | Distribution of forces between synergistics and antagonistics muscles using an optimization criterion depending on muscle contraction behavior[END_REF]. The use of a force-length relationship does not change the MusIC method structure. However, the study here should be extended to take into account the influence of the database density on the muscles length computation. The use of a force-velocity relationship would modify the MusIC method since the joint configuration space and the joint velocity space would have to be discretized for each sub-database. Thus, the study here should be extended to the muscles velocity computation.

Last, an alternative database generation method might be applied to decrease the database computation time without affecting the accuracy of the MusIC method. Indeed the discretization of the joint space may be non regular in order to favor the database density around joint configurations impacting more critically the solution. For example, moment arm sign changes are critical and need a high density of the database whereas extreme joint positions seem less sensitive and databases may be less dense at these locations. Such non regular discretization may be generated a priori for a given generic model by analyzing the sensitivity of moment arms changes with regard to joint configuration changes. However, the current study showed that accurate results can be obtained from low database densities. This result may mitigate the impact of a non regular discretization for models of this complexity. Non regular discretization may only be relevant for much more complex models.

Conclusion

The estimation of muscle forces using a precomputed database demonstrated significant computational benefits [START_REF] Muller | The MusIC method: a fast and quasi-optimal solution to the muscle forces estimation problem[END_REF]) with a minimal loss of accuracy. The current paper aimed at addressing the limitations of the MusIC method, mostly concerning the off-line computation time, to be used as a feedback during a musculoskeletal analysis session. The study modulated the database density to find a tradeoff between the accuracy of the method and the off-line computation time (database generation time) and validated the method on a 3-dimensions model currently used in musculoskeletal analysis. The method was developed to find a priori -without generating databases -the best trade-offs to be applied for a given model and ex-perimental design. It consists in using a generic model to evaluate the error made on moment arms computation during the joint space discretization. On the lower limbs musculoskeletal model available in the AnyBody Modelling System [START_REF] Damsgaard | Analysis of musculoskeletal systems in the AnyBody Modeling System[END_REF], the optimal database density allows to obtain an accuracy error of 0.1% with a corresponding off-line computation time lower than 10 minutes on a standard laptop. The on-line computation frequency is about 58 Hz. Therefore, thanks to this complementary work, the MusIC method may be used to provide a musculoskeletal analysis feedback in a really short time during a mocap-based musculoskeletal analysis session. In conclusion, this kind of method could allow musculoskeletal indicators to be integrated to an exercise or work pattern for sports training or ergonomics applications respectively.

Figure 2 .

 2 Figure 2. Databases structure composed of n j independent sub-databases (left part). Each of them is discretized at different joint configurations (central part) where the solution of the optimization problem is stored (right part).

Figure 3 .

 3 Figure 3. Representation of the musculoskeletal model used in this application. The model exhibits 12 DOF and 82 muscles.

Figure 4 .

 4 Figure 4. Different postures of the video guide. The video was generated using an avatar animated from motion capture data.

Figure 5 .

 5 Figure 5. Moment arm error ∆R dis for each joint space discretization. To facilitate the graph interpretation,