Victor Pillac

Christelle Guéret
email: gueret@mines-nantes.fr

Andrés L Medaglia

An event-driven optimization framework for dynamic vehicle routing

Keywords: dynamic vehicle routing, event-driven framework, multiple scenario approach, online stochastic optimization, DVRPSD, DVRP

The real-time operation of a fleet of vehicles introduces challenging optimization problems researches in a wide range of applications, thus, it is appealing to both academia and practitioners in industry. In this work we focus on dynamic vehicle routing problems and present an event-driven framework that can anticipate unknown changes in the problem information. The proposed framework is intrinsically parallelized to take advantage of modern multi-core and multi-threaded computing architectures. It is also designed to be easily embeddable in decision support systems that cope with a wide range of contexts and side constraints. We illustrate the flexibility of the framework by showing how it can be adapted to tackle the dynamic vehicle routing problem with stochastic demands. Computational results show that while our approach is competitive against state-of-the art algorithms, it still ensures greater reactivity and requires less assumptions (e.g., demand distributions).

Introduction

The problem of operating a fleet of vehicles arises in many contexts, from pickup and delivery of goods to relocation of trucks in carrier companies. More specifically, Vehicle Routing Problems (VRPs) deal with the design of a set of minimal-cost vehicle routes that serve the demand for goods or services of a group of geographically spread customers, satisfying a group of operational constraints. From an information perspective, such problems generally include two dimensions: evolution and quality of information [START_REF] Psaraftis | A dynamic-programming solution to the single vehicle many-to-many immediate request dial-a-ride problem[END_REF]. Information evolution relates to the fact that in some problems the information available to the planner may change during the execution of the routes, for example, with the arrival of new customer requests. Information quality reflects possible uncertainty on the available data, for instance, when the demand of a customer is only known as a range estimate of its real demand. In addition, depending on the problem and the available technology, vehicle routes can either be designed a-priori or online. Based on these dimensions, Table 1 identifies four categories of routing problems.

Information quality

Deterministic input Stochastic input

Information evolution

Input known beforehand Static and deterministic Static and stochastic

Input changes over time Dynamic and deterministic Dynamic and stochastic

Table 1: Taxonomy of vehicle routing problems by information evolution and quality.

The static and deterministic category includes the classical Vehicle Routing Problem (VRP) as defined by Dantzig [START_REF] Dantzig | The truck dispatching problem[END_REF] in which all information is known beforehand and with certainty. In contrast, problems from the static and stochastic class are characterized by input partially known as random variables, which realizations are only revealed during the execution of the routes. Additionally, it is assumed that routes are designed a-priori and only minor changes are allowed afterwards. A common example is the VRP with Stochastic Demands (VRPSD), in which customer demands are uncertain. We refer the interested reader to the surveys by Cordeau et al. [START_REF] Cordeau | Vehicle routing[END_REF], Baldacci et al. [START_REF] Baldacci | Recent advances in vehicle routing exact algorithms[END_REF], Laporte [START_REF] Laporte | What you should know about the vehicle routing problem[END_REF][START_REF] Laporte | Fifty years of vehicle routing[END_REF], and Toth and Vigo [START_REF] Toth | The vehicle routing problem[END_REF], for a recent review of these two classes of problems.

In dynamic and deterministic problems, also referred to as online problems, part or all of the input is unknown and revealed dynamically and unpredictably during the design or execution of the routes. On the other hand, dynamic and stochastic problems include partial stochastic knowledge on the dynamically revealed information. For these problems, vehicle routes are redefined in an ongoing fashion, requiring technological support for real time communication between the vehicles and the decision maker (e.g., mobile phones and global positioning systems). Techniques for both classes are reviewed in the studies by Ichoua et al. [START_REF] Ichoua | Planned route optimization for real-time vehicle routing[END_REF] and Pillac et al. [START_REF] Pillac | Dynamic vehicle routing: State of the art and prospects[END_REF].

Dynamism in routing can emerge from different aspects of the problem. The most common source of dynamism is the arrival of new customers with a demand for goods [START_REF] Hvattum | A branch-and-regret heuristic for stochastic and dynamic vehicle routing problems[END_REF][START_REF] Ichoua | Exploiting knowledge about future demands for real-time vehicle dispatching[END_REF][START_REF] Mes | Comparison of agent-based scheduling to look-ahead heuristics for real-time transportation problems[END_REF] or services [START_REF] Beaudry | Dynamic transportation of patients in hospitals[END_REF][START_REF] Gendreau | Parallel tabu search for real-time vehicle routing and dispatching[END_REF][START_REF] Goel | A general vehicle routing problem[END_REF]. Other researchers consider dynamically revealed demands for a set of known customers [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF][START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF][START_REF] Secomandi | Reoptimization approaches for the vehicle-routing problem with stochastic demands[END_REF], dynamic travel times [START_REF] Lorini | Online vehicle routing and scheduling with dynamic travel times[END_REF][START_REF] Potvin | Vehicle routing and scheduling with dynamic travel times[END_REF][START_REF] Tagmouti | A dynamic capacitated arc routing problem with time-dependent service costs[END_REF], and vehicle availability [START_REF] Li | A decision support system for the singledepot vehicle rescheduling problem[END_REF][START_REF] Li | Real-time vehicle rerouting problems with time windows[END_REF].

Fig. 1 illustrates the Dynamic Vehicle Routing Problem (DVRP), in which new customers appear while the vehicle is executing its route. Before the vehicle leaves the depot (at time t 0 ), an initial route plans to visit the currently known customers (A, B, C, D, E). While the vehicle executes its route, two new customers (X and Y ) appear at time t 1 and the initial route is adjusted to accommodate them. Finally (at time t f ), the executed route is (A, B, C, D, Y, E, X). This example reveals that dynamic routing requires to adjust the routes in an ongoing fashion, which implies real-time communication between vehicles and the dispatching center.

Until recently, the lack or high costs of real-time communication technologies steered vehicle routing research away from dynamic problems [START_REF] Eksioglu | The vehicle routing problem: A taxonomic review[END_REF]. Nevertheless, recent advances in communication and geolocation technologies now allow companies to economically track their fleet in real time. These new technologies lead to the development of Intelligent Transport Systems (ITS), The development of ITS and AFMS creates new challenges and opportunities for operations research. The advent of these systems demands a new class of efficient optimization algorithms to handle various difficult aspects of fleet management. Nevertheless, Crainic et al. [START_REF] Crainic | Intelligent freight-transportation systems: Assessment and the contribution of operations research[END_REF] suggest that while the hardware part of ITS has considerably evolved, the corresponding Decision Support Systems (DSSs), and optimization models in particular, have not yet reached their maturity.

From a practical perspective, we can identify the following desirable characteristics of a dynamic routing DSS:

• Event-driven. Advances in technology also mean that organizations are now able to react quickly to changes in their environment. Having a DSS which is periodically updated implies longer reaction delays. Thus, a DSS should be driven by the same transactional events that keep the business operating (e.g., customer requests).

• Parallelized. As dynamic routing requires fast decisions, the underlying optimization algorithms should be parallelized [START_REF] Crainic | Intelligent freight-transportation systems: Assessment and the contribution of operations research[END_REF][START_REF] Ichoua | Planned route optimization for real-time vehicle routing[END_REF], taking advantage of the now ubiquitous parallel (and distributed) computing architectures able to perform several tasks in parallel.

• Flexible. The landscape of vehicle routing problem variants is vast [START_REF] Eksioglu | The vehicle routing problem: A taxonomic review[END_REF]. Thus, a DSS should be easily extensible to account for specific aspects of different applications in a continuously evolving environment.

In this paper, we propose an application-oriented optimization framework for dynamic and stochastic vehicle routing that is event-driven, parallelized and flexible. The rest of this document is organized as follows. Section 2 reviews the literature on dynamic routing optimization techniques and related decision support systems. Section 3 describes the proposed framework, Section 4 illustrates its application to the dynamic VRPSD, and Section 5 presents experimental results. Finally, Section 6 concludes this paper and discusses how the framework can be generalized and extended to other dynamic optimization settings.

Literature review

A growing body of research has been carried out on dynamic routing, leading to new optimization techniques and innovative DSSs. In this section we will review some of the most significant contributions in the dynamic routing field.

Dynamic routing

A wide range of techniques have been developed to address the dynamic nature of routing problems. Dynamic methods can be divided in two categories: non-anticipative, which only react to updates in the problem data; and anticipative, which take into account knowledge on the dynamically revealed information to anticipate the future. Non-anticipative methods 2.1 Dynamic routing are designed for dynamic and deterministic problems. They generally are a direct adaptation of static methods such as integer programming [START_REF] Krumke | Real-time dispatching of guided and unguided automobile service units with soft time windows[END_REF][START_REF] Yang | Real-time multivehicle truckload pickup and delivery problems[END_REF], large neighborhood search [START_REF] Goel | A general vehicle routing problem[END_REF], tabu search [START_REF] Beaudry | Dynamic transportation of patients in hospitals[END_REF][START_REF] Gendreau | Parallel tabu search for real-time vehicle routing and dispatching[END_REF][START_REF] Ichoua | Vehicle dispatching with time-dependent travel times[END_REF], genetic algorithms [START_REF] Benyahia | Decision support for vehicle dispatching using genetic programming[END_REF][START_REF] Haghani | A dynamic vehicle routing problem with time-dependent travel times[END_REF][START_REF] Van Hemert | Dynamic routing problems with fruitful regions: Models and evolutionary computation[END_REF], and ant colony optimization [START_REF] Gambardella | Ant colony optimization for vehicle routing in advanced logistics systems[END_REF][START_REF] Montemanni | Ant colony system for a dynamic vehicle routing problem[END_REF]. Conversely, anticipative methods often make better decisions by using stochastic information available in the form of probability distributions. Anticipative methods are further classified into one of two families: stochastic modeling or sampling.

Anticipative methods based on stochastic modeling accurately describe the problem's stochasticity. In an early work, Powell et al. [START_REF] Powell | A comparative review of alternative algorithms for the dynamic vehicle allocation problem[END_REF] formulated the DVRP as a Markov Decision Process (MDP). Nevertheless, the exponential growth of the state and action spaces causes traditional MDPs to stall. This problem has led to the development of Approximate Dynamic Programming (ADP). The main idea behind ADP is to decompose the time in decision epochs. At each decision epoch the goal is to minimize the current deterministic cost plus an approximation of the expected future cost. This technique has been successfully applied to different dynamic fleet management problems [START_REF] Godfrey | An adaptive dynamic programming algorithm for dynamic fleet management, I: Single period travel times[END_REF][START_REF] Powell | Fleet management[END_REF][START_REF] Simao | An approximate dynamic programming algorithm for large-scale fleet management: A case application[END_REF], freight transport problems [START_REF] Powell | Approximate dynamic programming: solving the curses of dimensionality[END_REF][START_REF] Powell | Stochastic programming in transportation and logistics[END_REF], and vehicle routing with stochastic demands [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF]. The strength of ADP is that it accurately encapsulates stochastic information in the model, but at the expense of a higher complexity and stronger assumptions on the probability distributions.

On the other hand, anticipative methods based on sampling are to some extent simpler, but require more effort to capture the problem's stochasticity. These methods sample the probability distributions to generate scenarios that are used to make decisions. Such approaches include the dynamic sample scenario hedge heuristic proposed by Hvattum et al. [START_REF] Hvattum | Solving a dynamic and stochastic vehicle routing problem with a sample scenario hedging heuristic[END_REF], the tabu search heuristics proposed by Ichoua et al. [START_REF] Ichoua | Exploiting knowledge about future demands for real-time vehicle dispatching[END_REF] and Attanasio et al. [START_REF] Attanasio | Real-time fleet management at Ecourier Ltd[END_REF], and the Multiple Scenario Approach (MSA) proposed by Bent and Van Hentenryck [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF].

Among the anticipative methods based on sampling, MSA is unique in the sense that it provides a more general framework for dynamic problems. In particular, MSA maintains a pool of scenarios with realizations of the problem random variables and a solution to the corresponding deterministic problem. A distinctive feature of MSA is that the next customer to visit is selected based on the whole scenario pool by means of a decision process. Fig. 2 The algorithm starts by initializing the scenario pool based on the currently known information. At time step t, MSA updates the scenario pool to reflect the current environment state, selects the next customer, and optimizes the scenarios. As new information is disclosed, some scenarios might become obsolete and are removed from the pool, leaving space to new ones. The strength of MSA is that optimization is performed on scenarios and only requires to solve a static and deterministic problem. Therefore this approach is very flexible as it can 2.2 Decision support systems for dynamic routing virtually be adapted to any problem, provided an optimization algorithm for its static and deterministic version. Nonetheless, its integration in a real-world context is far from trivial, especially considering communication between the method and its environment. Additionally, the fact that it relies on time epochs induces delays between the arrival of new information and its processing.

Decision support systems for dynamic routing

Zak [START_REF] Zak | Decision support systems in transportation[END_REF] surveyed a wide range of DSSs for static vehicle routing. With our focus being on dynamic routing, this section reviews the body of research in this area.

The operation of a fleet of vehicles in an urban area is a key component of city logistics [START_REF] Taniguchi | City Logistics: Network Modelling and Intelligent Transport Systems. Pergamon[END_REF], and the core subject of various DSSs developments. For instance, Fleischmann et al. [START_REF] Fleischmann | Dynamic vehicle routing based on online traffic information[END_REF] presented an event-based DSS that takes into account changing travel times and the arrival of new customers in the context of a local area courier service. The framework continuously optimizes a single routing plan in which new customers are inserted either with an assignment model or insertion algorithms. A similar problem was addressed by Attanasio et al. [START_REF] Attanasio | Real-time fleet management at Ecourier Ltd[END_REF] who showed that the proposed DSS allows for an efficient operation (low administrative cost) as the fleet size (couriers) increases, a key competitive advantage in this sector. Likewise, Barcelo et al. [START_REF] Barcelo | Vehicle routing and scheduling models, simulation and city logistics[END_REF] presented a flexible DSS for vehicle routing and scheduling in city logistics and its application to the delivery of goods in two Italian cities. Their DSS includes a real time traffic simulator, connection to common GIS systems, and various routing models and optimization modules. In a different context, Zeimpekis et al. [START_REF] Zeimpekis | Dynamic management of a delayed delivery vehicle in a city logistics environment[END_REF] developed a DSS that takes into account unexpected events such as traffic conditions or vehicle breakdowns to re-optimize an existing distribution schedule. Li et al. [START_REF] Li | A decision support system for the singledepot vehicle rescheduling problem[END_REF] also studied vehicle breakdowns in an application to waste collection in Brazil. Dynamic DSSs generally rely on specific technology to ensure the communication between vehicles and the dispatching center [START_REF] Zeimpekis | Dynamic management of a delayed delivery vehicle in a city logistics environment[END_REF]. In contrast Bieding et al. [START_REF] Bieding | On line routing per mobile phone : A case on subsequent deliveries of newspapers[END_REF] propose a DSS based on a WAP (Wireless Application Protocol) server and mobile phones to manage the delivery of newspapers.

As pointed out by Crainic et al. [START_REF] Crainic | Intelligent freight-transportation systems: Assessment and the contribution of operations research[END_REF], there is a gap between state-of-the-art optimization techniques and the optimizers embedded in real-life DSSs. This may be explained by the complexity and level of specialization of certain approaches, that render difficult their extension and integration in an application-oriented context. To address this issue, we propose a flexible optimization framework, based on MSA [START_REF] Van Hentenryck | Online stochastic combinatorial optimization[END_REF], easily embeddable in any DSS for dynamic routing.

Proposed framework

The framework, called jMSA, is a flexible, parallel, and event-driven Java implementation of the multiple scenario approach. The proposed framework has been designed to facilitate and accelerate the development and deployment of MSA-based algorithms embeddable in DSSs. This section presents the proposed framework in detail.

Scenarios and decisions

Scenarios capture uncertainty in MSA. Each scenario contains a realization of the random variables, and a solution to the static and deterministic problem defined by this realization. For instance, in the Dynamic VRPSD (DVRPSD), in which vehicles can be dynamically rerouted, each scenario contains a realization of the customer demands; while in the DVRP, it contains a set of sampled (potential) customers, aside from the known customers. An optimization algorithm is used to solve the static and deterministic routing problem defined by both actual and sampled data. Fig. 3 illustrates how scenarios are generated for the DVRP. Solely based on the actual customers, the optimal tour would be (A, B, E, D, C), which ignores two zones (gray areas) where customers are likely to appear. By sampling the customer spatial distributions, customers X, Y and Z are generated, and the new optimal tour is (C, X, Y, B, A, Z, E, D). Removing the sampled (potential) customers leads to the tour (C, B, A, E, D) which is suboptimal based on a myopic cost evaluation, but leaves room to accommodate new customers at a lower cost. Another key element in MSA is the decision process, which defines how to select the next customer to serve based on the information of the scenario pool. MSA's accuracy relies to a great extent on the decision process, being the most common algorithms expectation, consensus, and regret. The expectation algorithm [START_REF] Chang | Online scheduling via sampling[END_REF] evaluates the cost of visiting each customer first, by forcing its visit and reoptimizing each scenario. The consensus algorithm [START_REF] Bent | Scenario-based planning for partially dynamic vehicle routing with stochastic customers[END_REF] selects the customer appearing first with the highest frequency. Finally, the regret algorithm [START_REF] Bent | Regrets only! online stochastic optimization under time constraints[END_REF] approximates the cost of visiting each customer first.

The jMSA framework unifies these decision processes in the generic Algorithm 1, in which a subset of candidate customers (line 1) is evaluated against the scenario pool (line 6) to select the best one (line 9). The evaluation of each customer reflects how desirable it is to serve it first depending on the objective. In most routing problems, the customer with the highest evaluation should be the one that ensures the lowest expected routing distance when visited first.

Algorithm 1 A general algorithm for the decision process in jMSA Input: scenario pool P, set of pending customers R Output: r * the next customer to serve

1: C ← selectCandidates(R, P) 2: f * ← -∞, r * ← ∅ 3: for all r ∈ C do 4: f ← 0 5:
for all s ∈ P do 6:

f ← f + evaluateRequestProfit(r, s) 7:
end for

8: if f > f * then 9: f * ← f, r * ← r 10:
end if 11: end for 12: return r *

Event-driven interaction

The original description of MSA is implicitly based on the notion of time steps introduced by Kilby et al. [START_REF] Kilby | Dynamic VRPs: a study of scenarios[END_REF]. By design, this time discretization in intervals implies a time lag between an update in the problem data, such as the arrival of a new customer, and the response of the system, corresponding to the time before the next time step. Consequently, in jMSA we propose a description of MSA from an event-driven perspective, suitable for its integration as a component of a real-world decision support system.

Fig. 4 illustrates a typical sequence of events while routing a single vehicle in a dynamic context. The environment refers to the real-world, the DSS is assumed to be based on the MSA algorithm, and active (idle) times are represented with a continuous (dotted) segment. While the vehicle is parked at the depot, the MSA procedure initializes a scenario pool based on the When the vehicle reaches its destination, an event is sent to the system (second dotted arrow) and triggers an update of the scenario pool. The remaining service time is used by MSA to reoptimize the pool until the vehicle is ready to depart. This event (third dotted arrow) triggers the decision procedure, which recommends visiting customer B (second double-headed arrow). At some point in time while the vehicle is traveling to B, an event (last dotted arrow) triggers an update of the scenario pool. Such event could be the arrival of a new customer in the DVRP, or an update in the traffic information in the case of routing with dynamic travel times.

The main advantage of this event-driven interaction between the environment and the system is that it increases the responsiveness of the DSS by feeding real-time information to the system and communicating decisions without delay.

Framework design

As illustrated in Fig. 5, the proposed framework is divided in two layers: a kernel, common to all dynamic combinatorial optimization problems; and a problem layer, with problem-specific components. The central component of the kernel is the MSAProcedure, which contains the logic of the algorithm and instantiates all other components. The MSAProcedure is configured via the Global-Parameters that can be set programmatically or via a configuration file.

Problem layer

Kernel

MSA Procedure

Event

The event-driven behavior is modeled using two elements: events and event handlers. Fig. 6 shows how events drive the framework. The MSA procedure continuously dequeues events from the event queue, and then processes them by using the corresponding event handler in the event handler manager.

Events are designed to increase the framework responsiveness. To ensure that important events are handled first, events are prioritized and the event queue is sorted accordingly. Additionally, some events are preemptive, meaning that handling of a non-preemptive event is always aborted in favor of a preemptive event. Event handlers define at a very high level what actions are triggered by a given event. By design, these handlers do not contain any problem-specific logic which is rather delegated to components. The component manager contains references to all components and acts as an interface between event handlers and problem-specific implementations. Fig. 7 illustrates how event handlers and components interact by means of the Scenario-Generation event. First, GenerateHandler calls the generateScenario method of the Component-Manager that internally uses the registered ScenarioGenerator. Then it calls the optimize-Scenario method, delegated to the instance of ScenarioOptimizer in use, and adds the scenario to the pool. The process repeats until the pool is full, moment when the event handling terminates by raising a ScenarioOptimization event that is further pushed to the event queue.

The framework includes a callback system that provides users with further control over the MSA procedure. Users may implement a callback simply by extending the generic Callback provided in the framework, and registering it in the MSA procedure. User-defined callbacks are automatically invoked at specific points of the procedure and allow customized uses such as logging to a file or dynamic parameter tunning.

Tied, yet decoupled to the kernel, the jMSA framework offers a problem-specific layer containing components that provide ready-to-use functionalities for common dynamic combinatorial optimization problems. Fig. 5 illustrates some components that could be combined for the DVRP.

Consensus is an implementation of the consensus algorithm that is common to many dynamic problems solved under MSA; VRPScenario is an implementation of Scenario for routing problems containing a set of routes; VRPScenarioOptimizer is a generic solver for the VRP; and finally DVRPScenarioGenerator is the only component specific to the DVRP that is responsible for the generation of new scenarios. This two-layer architecture ensures flexibility and extensibility. While kernel elements are defined at a high level and are designed to be problem independent, the problem layer provides implementations for specific problems. Thus, users only have to define or extend components, in particular for scenario generation and optimization, without worrying how they will be integrated in the MSA procedure.

Parallelization via multi-threading

The ubiquitous presence of multi-core processors can be exploited in parallelizable algorithms such as MSA. Nevertheless, parallelization often comes at the price of a higher (implementation) complexity. The jMSA framework offers multi-threaded parallelization of the most timeconsuming tasks, hiding it from the user. That is, under jMSA, users do not have to explicitly write a parallel algorithm, but simply rely on the ComponentManager which internally distributes tasks among different threads.

Fig. 8 illustrates how threads interact within the jMSA framework. At time t 0 the MSA thread dequeues an OptimizePool event, and processes it with the corresponding OptimizeHandler. In parallel to the MSA thread, two other threads are started by the ComponentManager to optimize the scenarios of the pool. At t 1 , a preemptive NewCustomer event is pushed by the environment, causing the MSA thread to prematurely abort the optimization. To avoid inconsistencies, the main thread waits for the pool executor to terminate, sends a signal to the callback thread to notify that the OptimizePool event was handled, and raises a GenerateScenarios event. Finally, the procedure dequeues the NewCustomer event, which has a higher priority than the Decision event, and processes it.

It is worth noting that aside from time-consuming tasks such as scenario generation and optimization, parallelization is also used to execute callbacks. Callbacks can be particularly useful when writing files or updating the state of a user interface as it does not affect the performance of the main algorithm. This behavior can be overridden using synchronous callbacks.

Problem description

The VRPSD consists in designing a-priori a minimal expected-cost routing plan composed of K routes, such that each customer is visited exactly once. The fundamental difference between the classic VRP and the VRPSD is that in the latter customer demands are known as random variables. The randomness in the VRPSD implies that a customer demand realization might exceed the vehicle remaining capacity, leading to a route failure that requires a recourse action.

An intuitive recourse action is for the vehicle to go back to the depot to restore its initial capacity and then resume its route [START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF], or to allow the service of additional customers before returning to the depot [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF]. It is important to stress that in this context all customers are known beforehand and the only dynamically revealed information is the realization of the customer demands. Uncertainty in the VRPSD has been addressed by various solution approaches, of which the two most studied are the Chance Constrained Programming (CCP) and the Stochastic Program-4.2 Scenarios and decisions ming with Recourse (SPR). Both methods are based on a two-stage approach: the first phase, builds robust routing plan; while the second phase takes place, takes recourse (corrective) actions as the realizations of the customer demands are unveiled. The conceptual difference between the two approaches lies in the objective of the first-stage optimization: in CCP, the goal is to ensure an upper bound on the probability of a failure, regardless of the expected cost of the second phase; while SPR seeks the minimization of the total expected cost, including recourse actions.

Henceforth, we focus on the single-vehicle Dynamic VRPSD (DVRPSD) (K = 1), where it is possible to reroute the vehicle upon new demand realizations, allowing more complex recourse actions. Literature on the DVRPSD is scarce, with the main contributions being the work by Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF], Secomandi [START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF], and Secomandi and Margot [START_REF] Secomandi | Reoptimization approaches for the vehicle-routing problem with stochastic demands[END_REF]. All of them consider the case where customer demands are discrete and uniformly distributed. We now show how the jMSA framework can be adapted to tackle this problem. Also, we illustrate how under the proposed approach we can easily relax the assumptions on the demand distributions, thus leading us to the solution of a more general problem with jMSA.

Scenarios and decisions

In the context of the DVRPSD, the only unknown data is the customer demand realization. Thus, scenarios contain different realizations of the customer demands, along with a feasible routing for these values. As the vehicle can go back to the depot during its service, a scenario can contain different routes that will be executed in a sequential order by the same vehicle.

The fact that customer locations are identical across scenarios suggests that different scenarios might have similar routes. Thus, we decided to use the consensus algorithm to select the next customer to visit. Let us consider the scenario pool of Fig. 9. The customers who have already been served (4 and 1) appear first in all scenarios, while customers 2, 3, 5, and 6, appear in varying order depending on the scenario sampled demands. Considering that customer 2 appears first in 2 out of 4 scenarios, by consensus it is selected as the next customer to visit. As shown in Algorithm 1, the function selectCandidates (line 1) returns the set of unserved customers while evaluateRequestProfit (line 6) returns 1 if customer r appears first in the scenario; 0, otherwise.

Optimization

To optimize scenarios we use an Adaptive Variable Neighborhood Search (AVNS), which is an extension of the Variable Neighborhood Search (VNS) [START_REF] Mladenovic | Variable neighborhood search[END_REF]. The main difference between AVNS and VNS is that neighborhoods are not explored sequentially, but randomly selected weighted by their previous performance. Our implementation uses an average ratio of the improvement to time as a metric of neighborhood performance, and maintains this information between calls of the optimization procedure. Neighborhoods with a better performance are more likely to be explored first, leading to a self-tuning algorithm. Our MSA scheme benefits from this automatic self-tuning behavior as the optimization procedure is called numerous times on similar instances (i.e., scenarios).

Algorithm 2 presents an outline of the AVNS algorithm. The algorithm initializes with the whole set of neighborhood structures (line 2), then it selects a neighborhood (line 4) to randomly perturb the current solution (line 5), and improves it by applying a local search procedure (line 6). If the new solution is accepted (line 8), then the new solution becomes the current solution N ← selectNeighborhood(N c ) 5:

x ′ ← shake(N, x)

6:

x ′ ← localSearch(x ′ )

7:

updatePerformance(N , x, x ′ ) 8:

if accept(x ′ , x) then 9:

x ← x ′ 10:

N c ← N 11:
else 12:

N c ← N c \ {N } 13:
end if 14:

if z(x ′ ) < z(x * ) then 15:

x * ← x ′ 16:

end if 17: end while 18: return x *

In our experiments we used the two neighborhoods structures Or-opt and string-exchange for the perturbation, and a Variable Neighborhood Descent (VND) based on swap and 2-opt as local search (line 6). A more detailed description of these neighborhoods can be found in the paper by Irnich et al. [START_REF] Irnich | Sequential search and its application to vehiclerouting problems[END_REF]. The initial solution is obtained by a Clarke and Wright (CW) heuristic [START_REF] Clarke | Scheduling of vehicles from a central depot to a number of delivery points[END_REF] in which the saving list is randomized, as presented in Mendoza et al. [START_REF] Mendoza | An evolutionary-based decision support system for vehicle routing: The case of a public utility[END_REF], leading to the CW+AVNS algorithm.

Failure handling

A route fails when a customer demand exceeds the vehicle's remaining capacity. Thus, the MSA procedure becomes aware of a route failure as soon as a Resource event is raised upon the arrival at the customer location. As a consequence, the route failures handling must be defined at the event handler level, by checking if the demand of the current customer is larger than the vehicle remaining capacity, and updating the scenario pool accordingly.

User interface

To illustrate the use of callbacks we developed a user interface shown in Fig. 10. The main panel (right) presents in real time the unserved (white) and served (dark gray) customers, the vehicle destination (light gray), and the executed route (arrows). The left panel, displays a log of events of jMSA and echoes the configuration settings. By means of a callback registered in the MSA procedure, all the information in the interface is updated in real time. The benchmark instances for the DVRPSD used in this work were initially proposed by Novoa [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF] and later used in Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF]. In this work we consider the larger problems with 30, 40 and 60 customers uniformly distributed in a 1 × 1 square grid with discrete uniform demands.

For each problem size, there are ten combinations of five different demand distributions by two vehicle capacities, leading to a complete testbed of 30 instances.

For each original instance (combination of problem size, demand distribution, and vehicle capacity) we generated 100 possible realizations of the customer demands, creating a complete testbed of 3000 instances. Each instance was solved to optimality using the COIN-OR Symphony VRP solver [START_REF] Ralphs | Symphony user manual[END_REF][START_REF] Ralphs | On the capacitated vehicle routing problem[END_REF].

CW+AVNS

To assess the optimization component in isolation we conducted an experiment on the 3000 instances derived from the Novoa [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF] benchmark. Fig. 11 presents the distribution of gaps to optimal values for the CW+AVNS algorithm and a CW+2-opt heuristic used as comparison. Note that CW+AVNS clearly dominates CW+2-opt, with 90% of all instances solved with a gap of less than 4%. Additionally, CW+AVNS runs relatively fast, with average CPU times between 50 ms and 650 ms for the larger instances.

DVRPSD

Our experimental setting is comparable to the one presented in Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF]. For each instance we ran a simulations using the jMSA framework as a black box. This means that an external simulator was used to send events to the MSA procedure simulating the vehicle route execution.

For easier comparison we report results in terms of value of information [START_REF] Mitrović-Minić | Waiting strategies for the dynamic pickup and delivery problem with time windows[END_REF]. The value of information for instance I, namely V(I), is the gap between the cost of the final solution returned [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF] instances by the algorithm z(I) and the a-posteriori optimal solution z * (I), and it is calculated as follows: Table 2 presents the results for problems from the Novoa testbed [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF]. Because of the aggregated results reported by Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF], we report the average value of information by using average solution values in Eq. 1. MSA dominates the algorithm proposed by Secomandi [START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF] (1s n2 r), and outperforms the best performing algorithms reported by Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF] (1s stostat r, 2s stostat r) for instances with 30 and 60 customers, and a vehicle capacity of 137 and 175. Additionally, MSA shows better overall results with an average gap of 3.3% against 4.8% for 2s stostat r, 5.8% for 1s stostat, and 13.6% for 1s n2 r. Aside from the performance in terms of value of information, it is important to stress that MSA runs continuously, and the next customer to visit is selected in a fraction of a second, while the other algorithms can take up to several minutes to make such decision, limiting their deployment and applicability in a real-world online DSS.

V(I) = z(I) -z * (I) z * (I) (1)
Aside from direct numerical comparison, the strength of our approach relies on the lack of strong assumptions on demand distributions. To illustrate this point, we adapted the testbed instances by converting demand discrete uniform distributions into left-truncated normal distributions (N LT ≥0 ) as follows:

U int (a, b) → N LT ≥0 a + b 2 , b -a + 2 6 (2)
Note that Eq. 2 ensures that the demand will be between a -1 and b + 1 with probability 0.997, and truncates negative values. as expected slightly better, with a reduction of 0.3% in the overall average value of information, which is due to the smaller variance. It is important to stress that to conduct this experiment the only change required was to use a different random number generator, which illustrates the flexibility of our approach. Other approaches based on stochastic modeling, like those of Novoa and Storer [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF], are not as flexible and heavily depend on distributional assumptions.

Conclusions

In this paper we presented the design and implementation of jMSA, an object-oriented eventdriven framework for the Multiple Scenario Approach (MSA). By doing a high-level abstraction of MSA to a problem independent level, we modeled it as an event-driven process that allows high reactivity to changes occurring in online and highly dynamic operational environments. We implemented jMSA as a flexible framework that is easily embeddable in decision support systems. By design, jMSA includes a callback system that gives the user further control over MSA and allow complex interactions with third party components. Additionally, we integrated into the framework the parallelization of time consuming tasks with no compromise for the framework user, which is a key aspect considering the wide availability of multi-core personal computers. We illustrated the use of jMSA on the DVRPSD. The optimization of scenarios is performed by an Adaptive Variable Neighborhood Search (AVNS) which improves an initial solution generated with a randomized Clarke and Wright heuristic. The strength of AVNS is that it automatically adjusts its search scheme depending on the problem's structure by keeping track of the neighborhood performance throughout the execution of the MSA procedure. Computational experiments show that our approach is competitive with state-of-the-art algorithms that take full advantage of the stochastic aspects, while it provides a more flexible scheme that can be used to tackle problems with more general demand distributions.

Figure 1 :

 1 Figure 1: Example of dynamic vehicle routing

Figure 2 :

 2 Figure 2: Overview of the original MSA algorithm flow diagram of MSA.The algorithm starts by initializing the scenario pool based on the currently known information. At time step t, MSA updates the scenario pool to reflect the current environment state, selects the next customer, and optimizes the scenarios. As new information is disclosed, some scenarios might become obsolete and are removed from the pool, leaving space to new ones.The strength of MSA is that optimization is performed on scenarios and only requires to solve a static and deterministic problem. Therefore this approach is very flexible as it can

Figure 3 :

 3 Figure 3: Scenario generation in MSA

Figure 4 :

 4 Figure 4: Time line of events for a the dynamic routing of a single vehicle

Figure 5 :

 5 Figure 5: Design overview of the jMSA framework

Figure 6 :

 6 Figure 6: Event-driven MSA framework

Figure 7 :

 7 Figure 7: Interaction between the GenerateHandler and the different components.

Figure 8 :

 8 Figure 8: Multiple threads interacting in jMSA

Figure 9 :

 9 Figure9: Example of the decision process by consensus in a 4-scenario pool. Each scenario contains customers who have been visited (in white) and customers yet to be visited (in gray).

Figure 10 :

 10 Figure 10: A graphical user interface for jMSA

  presents a high level

		Start
		(t=0)
	Initialize pool
	Update pool
		Decision	Next customer to visit
	Generate & optimize
		scenarios
		t=t+1
	NO	t=t f ?
		YES
		End

  , and the set of active neighborhood structures is reset (line 10). Otherwise, the current neighborhood is removed from the set of active neighborhoods (line 12). At each iteration, the performance of the current neighborhood is updated (line 7). This process iterates until all neighborhoods have been explored with no improvement.

			4.4 Failure handling
	(line 9)Algorithm 2 The Adaptive Variable Neighborhood Search algorithm
	Input: feasible solution x, evaluation function z , and set of neighborhood structures N =
	{N 1 , .., N K }
	Output: best solution found x *
	1: x * ← x
	2: N c ← N
	3: while N c = ∅ do
	4:
	École des Mines de Nantes	12/19	Technical Report 11/2/AUTO

Table 2 :

 2 Comparison of average value of information.

  Table 3 highlights the robustness of MSA which shows consistent performance when demand distributions are changed from uniform (discrete) to normal (continuous). Further, the results are

			Instance set (size,capacity)
	Algorithm	(30,137)	(30,87)	(40,183)	(40,116)	(60,274)	(60,175)	Average
	Uniform	0.9%	3.9%	3.5%	6.3%	2.9%	2.0%	3.3%
	Normal	0.7%	3.6%	3.4%	6.2%	2.2%	1.9%	3.0%

Table 3 :

 3 Comparison of average VI for discrete uniform and normal distributions.

École des Mines de Nantes 6/19 Technical Report 11/2/AUTO

Technical Report 11/2/AUTO

École des Mines de Nantes

École des Mines de Nantes 14/19 Technical Report 11/2/AUTO

École des Mines de Nantes 15/19 Technical Report 11/2/AUTO

École des Mines de Nantes 17/19 Technical Report 11/2/AUTO

Acknowledgements

Financial support for this work was provided by the Region Pays de la Loire (France), as part of the Vallée du Libre project, and the Centro de Estudios Interdisciplinarios Básicos y Aplicados en Complejidad (CeiBA, Colombia). This support is gratefully acknowledged.


Victor Pillac
email: vpillac@mines-nantes.fr

Christelle Guéret
email: gueret@mines-nantes.fr

Andrés L Medaglia
email: amedagli@uniandes.edu.co

A dynamic approach for the vehicle routing problem with stochastic demands

Keywords: dynamic vehicle routing, VRP with stochastic demands, multiple scenario approach

Introduction

The Vehicle Routing Problem with Stochastic Demands (VRPSD) is a variation of the classical Capacitated Vehicle Routing Problem (CVRP) in which the demand of each customer is modeled as a random variable and its realization is only known upon vehicle arrival to the customer site. Under this uncertain scenario, a possible outcome is that the demand of a customer ends up exceeding the remaining capacity of the vehicle, leading to a route failure. In this study we will focus on the single vehicle VRPSD in which the fleet is limited to one vehicle with finite capacity, that can execute various routes sequentially.

Traditional approaches for this problem involve a two-phase procedure [START_REF] Cordeau | Vehicle routing[END_REF][START_REF] Mendoza | A memetic algorithm for the multi-compartment vehicle routing problem with stochastic demands[END_REF]. In the first phase, the problem is solved a-priori with the goal of finding a set of robust routes ; while in the second phase routes failures are handled by means of recourse actions, such as forcing the vehicle to go back to the depot and restore its capacity, to then visit the client again and serve the remaining demand.

Recently, new approaches were developped based on the dynamic re-optimization of the vehicle routes [START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF][START_REF] Secomandi | Reoptimization approaches for the vehicle-routing problem with stochastic demands[END_REF][START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF]. In this context the routing is no longer done a-priori but instead performed (and re-optimized) throughout the day : as soon as the vehicle becomes idle, a decision has to be taken regarding the next customer visit. The algorithms proposed by the aforementioned authors offer greater flexibility in comparison with two-stage approaches, allowing significant cost reductions. However, they assume discrete uniform demand distributions and have to be run every time a decision has to be taken, which can be time consuming for large instances.

The present work is based on an adaptation of an optimization framework developed initially for the vehicle routing problem with dynamic customers (i.e., customers appear while the vehicles are executing their routes).

Approach

In the Multiple Scenario Approach (MSA) [START_REF] Van | Online stochastic combinatorial optimization[END_REF] a pool of scenarios is continuously optimized and updated to capture the uncertainties. Even if internally each scenario in the pool contains a complete routing for the customers, the only information visible to the user is the next customer to visit. This study proposes an adaptation of this framework to the VRPSD, by maintaining a pool of scenarios containing possible realizations of the customer demands. Scenarios are optimized with and Adaptive Variable Neighborhood Search and the next client to visit is selected with a Consensus algorithm.

Preliminary Results

Preliminary computational experiments show that our method is competitive against the algorithms proposed in [START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF][START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF], even if the test instances are not favorable as the customer locations are randomly distributed and their demands follow a uniform discrete distribution causing high variability in scenarios and realizations. TAB. 1 -Comparison of average gaps relative to the perfect information solutions Table 1 presents preliminary results for instances with 30, 40 and 60 customers from the Novoa benchmarks [START_REF] Novoa | Static and dynamic approaches for solving the vehicle routing problem with stochastic demands[END_REF] with an experimental setting identical to the one presented in [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF]. It can be noted that MSA dominates the algorithm proposed by Secomandi [START_REF] Secomandi | A rollout policy for the vehicle routing problem with stochastic demands[END_REF], while it outperforms the best performing algorithm reported in [START_REF] Novoa | An approximate dynamic programming approach for the vehicle routing problem with stochastic demands[END_REF] for the set of 60 customers instances with a vehicle capacity of 175. MSA is run continuously, and the next customer to visit is selected in a fraction of second, while the other algorithms can take up to various minutes.

Conclusions and research perspectives

We have demonstrated the flexibility of the Multiple Scenario Approach by adapting it to the VRPSD. Preliminary experiments show encouraging results, both in terms of routing cost and decision time, while relaxing assumptions on customer demand distributions. Further work will focus on the generalization to the multiple vehicles case and other customer demand distributions.


Victor Pillac
email: vpillac@mines-nantes.fr

Christelle Guéret
email: gueret@mines-nantes.fr

Andrés L Medaglia
email: amedagli@uniandes.edu.co

On the Technician Routing and Scheduling Problem

The technician routing and scheduling problem consists in routing and scheduling a crew of technicians in order to attend a set of service requests, subject to skill, tool, and spare part constraints.

In this study we propose a formal definition of the problem and present a constructive heuristic and a large neighborhood search optimization algorithm.

Introduction

The present work deals with the Technician Routing and Scheduling Problem (TRSP), in which a limited crew of technicians serves a set of requests. The TRSP can be seen as an extension of the Vehicle Routing Problem with Time Windows (VRPTW), where technicians play the role of vehicles and requests are made by clients. In the TRSP, each technician has a set of skills, tools, and spare parts, while requests require a subset of each. The problem is then to design a set of tours such that each request is visited exactly once, within its time window, by a technicians with the required skills, tools, and spare parts. The TRSP naturally arises in a wide range of applications, including telecoms, public utilities, and maintenance operations.

The TRSP is related to the multiple Traveling Salesman Problem with Time Windows (mTSPTW). At one end, if each technician has the required skills, tools, and spare parts to serve all requests, then the TRSP reduces to the mTSPTW. In contrast, if each request can be served by exactly one technician, the TRSP can be decomposed in independent TSPTWs (one per technician). The TRSP falls between the two problems, where both routing and assignment decisions are made, considering tools, spare parts, and return trips to the central depot.

The novelty of this problem is twofold. First, it introduces interesting compatibility constraints between technicians and requests. While skills are intrinsic attributes, technicians may carry different tools and spare parts over the planning horizon. Technicians usually start their tour from their home with an initial set of tools and spare parts that allows them to serve an initial set of requests. They also have the opportunity to replenish their tools and spare parts at a central depot at any time to service more requests. Tools can be seen as renewable resources, while spare parts are non renewable and consumed once the technician serves a customer. Second, the TRSP differs from most routing problems (including the classical VRP) in its objective, as the workload balance among technicians is at least as important as the minimization of the traveled distance.

In this work we focus on a slightly restricted version of the TRSP, where requests require one skill and have hard time windows. We also assume that while tools and spare parts are available in unlimited quantities at the central depot, technicians are indeed limited. The objective is to balance the workload, while meeting a maximum shift length.

Literature review

Despite its numerous applications, technician routing problems have received limited attention and relatively few research has been conducted on closely related problems.

Borenstein et al. [START_REF] Borenstein | On the partitioning of dynamic workforce scheduling problems[END_REF] studied the technician workforce scheduling problem faced by British Telecom considering dynamic requests and skill constraints. To reduce its complexity they group requests into clusters and heuristically assign technicians to areas considering multiple objectives. They propose a rule based system that assigns and schedules requests to technicians depending on metrics such as their priority.

Xu and Chiu [START_REF] Xu | Effective heuristic procedures for a field technician scheduling problem[END_REF] studied a variant of the TRSP in which the objective is to maximize a weighted sum of job priorities and the technicians' idle time. Tang et al. [START_REF] Tang | Scheduling technicians for planned maintenance of geographically distributed equipment[END_REF] also considered job priorities and formulated the TRSP as a multi-period maximum collection problem, in which time-dependent rewards are granted for the completion of a request.

Finally, Bredström and Rönnqvist [START_REF] Bredström | Combined vehicle routing and scheduling with temporal precedence and synchronization constraints[END_REF] studied an interesting variant of the VRPTW in which some clients must be visited simultaneously by two or more vehicles. Although they do not explicitly consider skills, the proposed model accounts for compatibility constraints between vehicles and requests.

To the best of our knowledge, no author has considered tools or spare parts. In the following section we will describe an approach to tackle this problem.

Proposed approach

The proposed approach couples two components, namely, a constructive heuristic with an Adaptive Large Neighborhood Search (ALNS) procedure.

The constructive heuristic is a parallelized greedy algorithm that constructs a tour for each technician by iteratively inserting requests. This is achieved by evaluating the cost of inserting all requests in all tours, and then solving an assignment problem to decide which request will be inserted in each tour. A direct consequence is that at each iteration the algorithm inserts at most one request in each tour. The procedure only considers feasible insertions, possibly inserting a trip to the depot to satisfy tool and spare part constraints. If no further requests can be inserted, the algorithm backtracks and reverts the first assignments, storing removed requests in a temporary pool to prevent cycling. Then it attempts the insertion of unserved requests. If successful, the procedure is resumed with the whole set of unserved requests; if not, additional requests are removed. The algorithm stops when a feasible solution is found or after a predefined number of iterations. If no feasible solution can be found, then the solution servicing the largest number of requests is used.

The ALNS algorithm [START_REF] Pisinger | Large neighborhood search[END_REF] works by successively destroying and repairing a current solution. At each iteration, the algorithm selects destroy and repair procedure with a roulette wheel that reflects their past performance. Destroy methods remove a subset of 10% to 40% of the requests from the current solution, while repair procedures reinsert them using heuristics that are known to perform well on the problem at hand. The resulting new solution is accepted as current solution according to a simulated annealing criterion. The algorithm stops after a fixed number of iterations.

We developed three destroy procedures for the ALNS: the first, random destroy, consists in randomly removing requests, allowing diversification of the search; the second, related destroy, selects requests that are related in terms of required skills, time windows, and location, the goal being to exchange similar requests among technicians; finally, critical destroy, attempts to removes requests that are costly in the current solution. In all three cases, the goal is to reverse assignment decisions that were taken previously and overconstrain the technician tours (e.g., by forcing the visit of requests in distant areas or with conflicting time windows). In addition, destroy procedures evaluate the removal of unnecessary trips to the central depot.

On the other hand, we use two repair procedures. The first one, best insertion, is an intuitive greedy heuristic that iteratively inserts the request with the lowest insertion cost. The second procedure, regret, estimates the cost of not inserting a request in its best tour, in order to avoid postponing difficult insertions, as it can occur with best insertion.

In addition, we propose a Variable Neighborhood Descent (VND) [START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF] that optimizes each tour independently, without considering exchanges of requests between technicians. Two neighborhoods are used: shift, which moves a request within the tour, and 2-opt; which reverses a subtour. At this stage, Udine, Italy, July 25-28, 2011 preventive trips to the central depot that were previously planned to ensure tour feasibility can be moved within the tour, provided that tool and spare parts constraints are satisfied. This cluster-first, route-second approach is justified by the strong compatibility constraints between requests and technicians that limit the possible intertour optimization, and by the fact that the destroy-repair process addresses this issue. The VND is used as a post-optimization procedure in both the constructive and optimization phases.

To assess the performance of the proposed approach we adapted the instances of Solomon [START_REF] Solomon | Algorithms for the vehicle-routing and scheduling problems with time window constraints[END_REF] adding skills, tools, and spare parts information. Preliminary results show that a basic version of the ALNS is able to improve the initial solution total working time by 5% and the standard deviation of tour length by 7%. In addition, work currently underway involves testing the algorithm's adaptability to problems faced by maintenance companies.

Conclusions

We described the TRSP, a routing problem with both challenging features and numerous applications. As a first step towards the dynamic version, this work proposes a metaheuristic approach that tackles the static version of the problem, in which all information is known beforehand. Future work will focus on the development of online algorithms able to cope with both planned (e.g., preventive maintenance) and dynamic requests (e.g., corrective maintenance), that may emerge during the day, as well as unexpected delays such as longer service times.

MIC 2011 :

 2011 The IX Metaheuristics International Conference S2-40-3

Acknowledgements

Financial support for this work was provided by the Region Pays de la Loire (France), as part of the Vallée du Libre project, and the Centro de Estudios Interdisciplinarios Básicos y Aplicados en Complejidad (CEIBA, Colombia). This support is gratefully acknowledged.


Sébastien Baguet
email: baguet@mn.esm2.imt-mrs.fr

Bruno Cochelin
email: cochelin@mn.esm2.imt-mrs.fr

Stability of thin-shell structures and imperfection sensitivity analysis with the Asymptotic Numerical Method

Keywords: Buckling, thin shells, imperfection sensitivity, fold line, extended system, Asymptotic Numerical Method, finite elements. MOTS-CLÉS : Flambement, coques minces, sensibilité aux imperfections, courbe de suivi de points limites, système augmenté, Méthode Asymptotique Numérique, éléments finis

This paper is concerned with stability behaviour and imperfection sensitivity of thin elastic shells. The aim is to determine the reduction of the critical buckling load as a function of the imperfection amplitude. For this purpose, the direct calculation of the so-called fold line connecting all the limit points of the equilibrium branches when the imperfection varies is performed. This fold line is the solution of an extended system demanding the criticality of the equilibrium. The Asymptotic Numerical Method is used as an alternative to Newton-like incremental-iterative procedures for solving this extended system. It results in a very robust and efficient path-following algorithm that takes the singularity of the tangent stiffness matrix into account. Two specific types of imperfections are detailed and several numerical examples are discussed.

RÉSUMÉ. Ce travail concerne la stabilité des coques minces élastiques et leur sensibilité aux imperfections. L'objectif est la détermination de la courbe de réduction de charge critique de flambement en fonction de l'amplitude de l'imperfection. Pour cela, on procède à un calcul direct de la courbe de suivi de point limites qui relie tous les points limites des branches d'équilibre lorsque l'imperfection varie. Cette courbe de suivi est solution d'un système augmenté qui caractérise les points d'équilibre singuliers. La résolution de ce système augmenté est effectuée à l'aide de la Méthode Asymptotique Numérique qui est une alternative aux méthodes incrémentales itératives de type Newton. L'algorithme de suivi de branche qui en résulte est très robuste et très performant et prend en compte la singularité de la matrice de rigidité tangente. Deux types d'imperfections sont détaillés et plusieurs exemples sont traités.

Introduction

The development of new materials and the wide use of optimization techniques over the past decades have conducted to very thin and slender structures. Since structural stability is a question of prime interest for such structures, a considerable attention has been given to the development of efficient numerical algorithms for the description of the post-buckling path and the treatment of bifurcation and limit points which provide the critical buckling load (see figure 1). In the case of quasi-static linear elasticity, the underlying post-buckling problem can be represented by the discretized geometrically non-linear equation

F (u, λ) = f (u) -p(λ) = 0 [1]
which expresses the balance between internal forces and external forces. The internal forces depend only on the current displacement vector u where u is a N -dimensional vector, i.e. contains N dof's. The external forces are assumed to be proportional to an external loading F e and take the form p(λ) = λF e where λ is the load parameter. The solution of [1] is usually obtained by an incremental-iterative procedure such as the Newton-Raphson method [RIK 72, CRI 81] which is commonly implemented in standard commercial FEM codes. However, alternative procedures for solving nonlinear problems in structural mechanics, like the Asymptotic Numerical Method (ANM) [DAM 90, AZR 93, COC 94b] or the LATIN method [LAD 98], can also be used.

For slender structures, a stability analysis is all the more complex that they are often very sensitive to imperfections. A small defect can significantly reduce their critical buckling load, causing the collapse for a smaller value than in the perfect case (without imperfection). A simple way for carrying a sensitivity analysis consists in introducing an imperfection in the definition of the structure and in solving problem [1] for this "imperfect" structure. Doing so, an equilibrium branch is obtained for each value of the imperfection amplitude Λ, as represented in figure 1, and the limit point on each of them gives the "imperfect" buckling load. This type of analysis is very simple to implement and the previously described algorithms can be used.

Nevertheless, this procedure is not the most suitable for a sensitivity analysis. Indeed, this succession of re-analyses is very costly since a full computation on the structure is required for each requested limit point. Moreover, it provides informations only for the few imperfection amplitudes that have been fixed. In practice, it is not judicious to calculate all the different equilibrium paths. It is more advisable to determine precisely a starting critical state for a given imperfection amplitude, and then to follow directly the fold line connecting all the limit points when the imperfection varies (see figure 1). By this way, it is not necessary to compute all the equilibrium paths. Only the fold line is computed. To this end, the imperfection amplitude Λ must become an additional parameter in the governing equations. In this paper we will only consider shape and thickness imperfections . Thus, only the internal forces are affected and the previous one-parameter system [1] becomes a two-parameter system which reads

F (u, λ, Λ) = f (u, Λ) -λF e = 0 [2]
The solution of this system is now a load-imperfection-displacement surface which is made of the load-displacement curves for all the possible values of the imperfection amplitude. On this surface, we are only interested in the limit points with respect to λ. In order to isolate these limit points, a constraint equation characterizing the criticality of the equilibrium is appended to the equilibrium equations. The two-parameter nonlinear system [2] is then transformed into a so-called "augmented" or "extended" system which reads

H (u, λ, Λ) = F (u, λ, Λ) G (u, Λ) = 0 [3]
The solution of this extended system is restricted to the fold-line that directly gives the reduction of the critical buckling load when the imperfection varies.

Augmented problem

Many alternatives concerning the constraint equation G (u, Λ) = 0 have been proposed in the literature. The simplest of them is based on the study of the determinant of the tangent stiffness matrix K T = F ,u , where the notation F ,u indicates the differentiation of F with respect to u, but this criterium is not well suited for a numerical implementation. Another criterium which is more efficient from a numerical point of view is also based on the appearance of a null eigenvalue for the tangent stiffness matrix K T at simple critical states but it does not require the computation of the determinant. The right null vector φ associated with the null eigenvalue of K T is used instead, yielding the following constraint equation

G (u, Λ) = F ,u (u, Λ, λ) • φ = K T (u, Λ) • φ = 0 [4]
With this expression for the additional constraint equation, the extended system [3] can be rewritten as

H (u, φ, Λ, λ) =   F (u, Λ, λ) F ,u (u, Λ, λ) • φ φ -1   =   f (u, Λ) -λF e K T (u, Λ) • φ φ -1   = 0 [5]
The last normalization equation on φ is added to ensure the uniqueness of the solution and to avoid numerical problems by preventing φ from becoming excessively large or small.

In the context of structural mechanics, this extended system has been first introduced by Keener

Asymptotic algorithm for the fold line following

In this paper, the sensitivity analysis is reconsidered using the so-called Asymptotic Numerical Method (ANM) as an alternative to classical incremental-iterative methods. This method is inspired by the perturbation techniques which were developed by Thompson and Walker [THO 68] for decomposing a nonlinear problem into a sequence of linear ones. These perturbation techniques have been revisited and efficiently solved by means of the finite element method by Cochelin, Damil and Potier-Ferry [DAM 90, AZR 93, COC 94b]. Since the solution obtained with the ANM is valid only in the vicinity of the starting point, a continuation method has been proposed by Cochelin [COC 94a], allowing a curve to be described in a step by step way.

The ANM has two major advantages. On one hand, its computational cost is reduced as compared to the incremental-iterative methods because the continuation procedure requires less tangent stiffness matrix decompositions. On the other hand, the continuation method is very robust and reliable. It allows to follow very complicated paths without jumping on other branches as it can happen with the Newton-Raphson method when the length of the predictor step is not carefully set. With the ANM, the step length is based on the equilibrium residual and can be computed a posteriori. By this way, it is guaranteed to be always optimal for each continuation step and no special procedure is needed for its resizing. This method can be used for solving a large class of problems involving smooth or strong nonlinearities. For a non-exhaustive list, the interested reader should refer to [POT 97,NAJ 98]. An original bifurcation indicator based on the ANM has been proposed by Boutyour [BOU 95, BAG 00] for the precise detection of simple bifurcations and limit points and for the calculation of the associated eigenmode. Supposing that a starting limit point has been isolated by means of this algorithm, the direct following of the fold line is performed as follows with the ANM.

We assume that the solution (u, φ, λ, Λ) of the extended system [5] can be represented by a truncated power series expansion with respect to an additional parameter a

u(a) = u 0 + a u 1 + a 2 u 2 + . . . + a n u n φ(a) = φ 0 + a φ 1 + a 2 φ 2 + . . . + a n φ n Λ(a) = Λ 0 + a Λ 1 + a 2 Λ 2 + . . . + a n Λ n λ(a) = λ 0 + a λ 1 + a 2 λ 2 + . . . + a n λ n [6]
where the point (u 0 , φ 0 , Λ 0 , λ 0 ) is known. In fact, it corresponds to the previously isolated limit point which is used as a starting point. This starting point is a singular point of the equilibrium equation F , like every point of the fold line, but it is a regular point of the extended system H. These series expansions depend on the path parameter a. Because of this new variable, an additional constraint equation is needed. By analogy with the classical arc-length method, we choose the following definition for a

a = (u -u 0 ).u 1 + (Λ -Λ 0 )Λ 1 + (λ -λ 0 )λ 1 [7]
The next step consists in introducing the series expansions [6] into the nonlinear extended system [5] and in the definition of a [7]. By identifying the terms with the same power of a, one obtains a succession of linear problems and the problem at order

1 reads     K T 0 F 1 -F e K φ K T F 2 0 0 φ T 0 0 0 u T 1 0 Λ 1 λ 1         u 1 φ 1 Λ 1 λ 1     =     0 0 0 1    
(N eq.) (N eq.) (1 eq.) (1 eq.) [8] The solution (u 1 , φ 1 , Λ 1 , λ 1 ) of this linear system is the tangential direction of the fold line at the starting point (u 0 , φ 0 , Λ 0 , λ 0 ). This system is exactly the same as for the predictor step of the Newton-Raphson method [ERI 99, DEM 97]. Indeed, the matrix K φ et the two vectors F 1 and F 2 correspond to the differentials of the mapping F evaluated at the starting point (u 0 , φ 0 , Λ 0 , λ 0 )

K T (u 0 , Λ 0 ) = f ,u (u 0 , Λ 0 ) [9] K φ (u 0 , φ 0 , Λ 0 ) = f ,uu (u 0 , Λ 0 ) • φ 0 [10] F 1 (u 0 , Λ 0 ) = f ,Λ (u 0 , Λ 0 ) [11] F 2 (u 0 , φ 0 , Λ 0 ) = f ,uΛ (u 0 , Λ 0 ) • φ 0 [12]
As pointed out by several authors [ERI 99], the accuracy of these derivatives are of paramount importance for the efficiency of the algorithm because they are used to define the properties of the solution. Nevertheless, since an analytical calculation of

Produit par HAL -16 Sep 2014

these differentials is often a very painful task, a numerical approximation is usually preferred [WRI 90, REI 95]. In the present case, with the use of the ANM, exact analytical forms of these differentials can be easily obtained. Since the introduction of series expansions into a set of equations acts as a linearization tool, the directional derivatives with respect to the expanded variables follow from the identification of the power-like terms. Moreover, provided that the equations are written in a quadratic framework, the asymptotic expansions can easily be computed. The different steps that lead to the extended system [8] at order 1 and provide K T , K φ , F 1 and F 2 are described in appendix I. It must also be noted that, because of the differentiation with respect to Λ, the vectors F 1 and F 2 depend on the type of imperfection (shape or thickness in the present paper) whereas the matrices K T and K φ do not.

Equating power-like terms at order p (p ≤ n) yields the following system

    K T 0 F 1 -F e K φ K T F 2 0 0 φ T 0 0 0 u T 1 0 Λ 1 λ 1         u p φ p Λ p λ p     =     F nl p G nl p h nl p 0    
(N eq.) (N eq.) (1 eq.) (1 eq.)

[13]

The tangent augmented matrix is the same as for order 1. Only the r.h.s. terms are different. These terms contain the non-linearities of the problem. They are fully determined because they depend only on the solution at previous orders (u r , φ r , λ r , Λ r ) with r < p. As a result, recursively solving the succession of linear augmented systems will provide all the coefficients of the series [6]. The r.h.s. terms depend on the considered type of imperfection and, as for the differentials [9]-[12], analytical expressions follow from the identification process. For lack of space, their expression will not be given here. The interested reader can refer to [BAG 01] for more details about their calculation.

In practice, a deflated block elimination is used to solve the extended system [13] in order to consider only subsystems of size N involving the matrix K T . Such a blockelimination scheme can be found in Wriggers and Simo [WRI 90]. Its main interest relies on the fact that only the classical matrix K T needs to be decomposed, thus saving a large amount of calculation time. Besides this particular procedure, another numerical difficulty must be pointed out. Since all the solution points of H are singular ones of F , the matrix K T is singular all along the fold line connecting the computed solution points. That means that the classical matrix decomposition techniques cannot be used. A special procedure based on Lagrange multipliers is introduced to bypass this problem. A detailed description of the adaptation of the block-elimination and Lagrange multipliers procedures to the system [13] is given in [BAG 00].

Because of the limited radius of convergence of the series, only a fraction of the solution curve is obtained. Indeed, the accuracy of the solution deteriorates very quickly when the radius of convergence is reached. A criterium that gives the length of a step has been proposed by Cochelin [COC 94a]. Its justification is given in [BAG 01]. This criterium is based on the study of the residual of equilibrium equations [2]. In the case of our extended system, both residuals of the first and second equations of [5] must be Produit par HAL -16 Sep 2014 monitored. However, numerical experiments have shown that the second equation is always more accurate than the first one. As a result, the classical criterium based on the equilibrium residual can be used. According to this criterium, for series truncated at order n, the maximal value of the path parameter a for which the solution satisfies a requested accuracy ε is given by

a M = ε F nl n+1 1 n+1 [14]
Using this formula, the step length is set after all the coefficients of the series have been computed. By this way, the step length is guaranteed to be optimal. Once the step has been stopped, the starting point is updated and the global procedure is restarted. Thus, the solution curve is described in a step by step way, as it would be with the classical continuation algorithms. The strong point of this procedure is its robustness. Furthermore, it is completely automatic from the user's point of view. The only parameters that need to be chosen are the order n of the series and the accuracy ε. Setting n equal to 20 or 30 and ε = 10 -5 is often a good compromise.

Numerical results

The following examples are concerned with a cylindrical panel submitted to a concentrated load at its top. The geometrical and material properties and the boundary conditions of the panel are defined in Figure 2a. Using symmetry conditions, only a quarter of the panel was considered for the calculations and a 10×10 mesh with 200 triangular DKT shell elements [BAT 92] and 726 degrees of freedom was used. The aim of these examples is to analyse the sensitivity of the critical buckling load of the cylindrical panel to a thickness or a shape imperfection, i.e. to calculate the fold line for both cases. Using the Asymptotic Numerical Method, the equilibrium branch was first evaluated for a thickness h=6.35mm. The well-known results of this academic test are reported in Figure 2b. The equilibrium branch exhibits two limit points (one maximum and one minimum) associated with the snap-through behaviour of the panel. These limit points can subsequently be used as starting points for the fold line following. With an accuracy ε = 10 -5 and series at order 30, this equilibrium curve was obtained with only 12 ANM continuation steps, i.e. with only 12 decompositions of the tangent stiffness matrix K T .

Cylindrical panel with thickness imperfection

In this first example, the thickness h is chosen as the additional parameter Λ in order to carry out a sensitivity analysis with respect to the thickness of the panel. This can be used for example to study the influence of corrosion on the stability of the structure. A starting limit point was first isolated for an initial value of the thickness h=12.7mm and then the fold line was followed using the Asymptotic Numerical Method as described in previous sections. The "load-deflection-thickness" 3D representation of this fold line as well as the associated projections are plotted in Figure 3. Only 20 ANM continuation steps with series at order 30, i.e. only 20 decompositions of the tangent stiffness matrix K T , were needed to describe the entire fold line. The end of each continuation step is represented by a dot in figure 3.

The two most interesting projections of the fold line are shown in figures 4 and 5. In figure 4 the fold line as well as two equilibrium curves (for h=12.7mm and h=22.15mm) with their associated limit points are plotted. As expected, the fold line intersects the equilibrium branches at their two limit points (maximum and minimum), excepted for the extreme value h≃22.15 of the thickness for which the two limit points merge into an inflection point and exchange their role. Beyond this value, there is no snap-through phenomenon anymore. This behaviour can also be observed on figure 5 which shows the critical buckling load as a function of the panel thickness. For a given thickness, this curve gives the critical buckling load corresponding to the maximum and the minimum limit points of the equilibrium branches of figure 4. It can also clearly be inferred from this curve that there is no limit point beyond the value h≃22.15. The results obtained here are in very good agreement with those obtained by Eriksson et al. in [ERI 99].

Most of the time, only a part of the fold line in the neighbourhood of a given imperfection amplitude is needed. In that case, it is not useful to compute the entire fold Produit par HAL -16 Sep 2014 line. The following strategy, as illustrated in figure 6, is preferable. A starting limit point was isolated for a value of the thickness close to h=22.15. Then the fold line was followed in the neighbourhood of this value. Since the series expansions [6] are valid for positive or negative values of the path parameter a, the fold line can be followed in the two directions with the same series. Moreover, the series expansions can beneficially be replaced by rational representations called Padé approximants. The extra calculation cost is very small and the range of validity of the series is usually doubled [COC 94c, ELH 00]. The results obtained with Padé approximants are plotted in figure [6]. The previously computed, complete fold line is also plotted as a reference. No continuation was used for Padé approximants, only one step is represented here, and no residual criterium was used to stop the step. So, one can see their range of validity. With only one decomposition of the tangent matrix K T and the use of Padé approximants, the fold line is valid from h≃22.15 to h≃12.7. Beyond this value, the solution is out of the range of validity and the curve deviates from the reference fold line.

Cylindrical panel with shape imperfection

The second example refers to the same cylindrical panel, with geometrical and material properties as defined in figure 2. The aim is now to analyse the sensitivity of its

Produit par HAL -16 Sep 2014

critical buckling load to a geometrical shape imperfection. The first buckling eigenmode of the panel is shown in Figure 7. It was normalized by setting the greatest value of its displacement to 1 in order to provide the shape of the considered geometrical imperfection. The additional scalar parameter Λ is now used to monitor the amplitude of the shape imperfection. The notation η will be used instead of Λ in the sequel. For this example, the starting point was evaluated for an initial imperfection amplitude η 0 =-3mm, i.e. -0.47 times the thickness (h=6.35mm) of the panel. The curve shown in figure 8 is a 3D representation of the resulting fold line connecting all the limit points when the imperfection varies. This is a closed curve. Its projections are more suitable for our analysis. Two of them are represented in figures 9 and 10. The first projection is located in the load-deflection plane. Three equilibrium branches are plotted for different values of the imperfection as well as the projection of the fold line. The fold line intersects each equilibrium branch at its two limit points. The equilibrium branch of figure 9 for the perfect panel (η=0) is represented here again as a reference. The values of the load at its two limit points (maximum and minimum) can be also read on the right-hand-side figure 10 which gives the critical buckling load reduction with respect to the imperfection. The top and bottom parts of this curve correspond respectively to the variation of the maximum and minimum limit point of the equilibrium curves. The curve shows that the two limit points get closer when the imperfection amplitude increases. For two extreme values η≃-3.5h and η≃5.3h of the imperfection, the limit points merge and exchange their role. Beyond these extreme values, there is no limit point anymore, i.e. no snap-through phenomenon anymore. It can be noted that the entire fold line was computed with only 25 continuation steps, i.e. only 25 tangent matrix decompositions.

Conclusions

This paper has described a general procedure for the imperfection sensitivity analysis of elastic structures. The developed calculation tool can handle global as well as local imperfections, in the case of thickness or geometrical shape defects.

The Asymptotic Numerical Method is used for the numerical treatment of the problem. As a result, the required calculation cost is significantly reduced as compared to the classical Newton-Raphson procedure. Moreover, the geometrical nonlinearities are treated without any approximation and exact analytical expressions are obtained

Produit par HAL -16 Sep 2014

for the directional derivatives of the tangent stiffness matrix. Thus, no approximation is introduced and the resulting overall algorithm is very accurate and the continuation procedure is very robust.

In order to deal with more complex structures, future developments will concern the implementation of a shell element which handles large rotations. Material nonlinearities will also be introduced. Nevertheless, because of the dependance to the loading history, plasticity with unloading phenomenon can not be considered for the fold line following and the only possible nonlinear constitutive laws are be restricted to nonlinear elasticity. In the case of plastic buckling, an alternative strategy to the fold line following must be used. Such a strategy can be found for example in [LEG 02].

The next step consists in introducing the series expansions [6] of u, S, φ, Ψ, η and λ into [18] and [19]. The identification of the terms with the same power of a gives at order 1

     Ω BT 0 S 1 + B T nl (u 1 ) S 0 + η 1 B nl (u ⋆ 0 ) S 0 dV = λ 1 F e S 1 = D B0 u 1 + η 1 B nl (u 0 ) u ⋆ 0 [20]      Ω BT 0 Ψ 1 + B T nl (φ 1 ) S 0 + B T nl (u 1 ) Ψ + η 1 B T nl (u ⋆ 0 ) Ψ 0 + B T nl (φ 0 ) S 1 dV = 0 Ψ 1 = D B0 φ 1 + B nl (φ 0 ) (u 1 + η 1 u ⋆ 0 )
[21] B0 = B(u 0 + η 0 u ⋆ 0 ) is introduced to shorten the notations. In order to go back to a displacement formulation, S 1 and Ψ 1 are then replaced in equations [20] 1 and [21] 1 . These variables were introduced only to make the asymptotic expansions easier. Finally, the problem at order 1 can be written in the following form

K T u p + η p F 1 = λ p F [22] K φ u p + K T φ p + η p F 2 = 0 [23]
which corresponds to the two first lines of [8] with In this section, the additional parameter is the thickness of the shell. For simplicity, we focus on shell finite elements which are analytically integrated through the thickness, and where the thickness explicitly appears only in the elasticity matrix D. The DKT18 shell element [BAT 92] is an example of such an element. For this class of shell elements, the elasticity matrix D(h) = D 1 h + D 3 h 3 can be decomposed into a membrane part and into a bending part, which are respectively linear and cubic with respect to h. According to [15] 2 and the definition of D, the constitutive law is quadratic with respect to u and cubic with respect to h. Since the asymptotic expansions are easier in a quadratic framework, we introduce some additional variables in Produit par HAL -16 Sep 2014

K t = Ω BT 0 D B0 + G T Ŝ0 G dV [24]

Figure 1 .

 1 Figure 1. Perfect and imperfect equilibrium paths and fold line

Figure 2 .

 2 Figure 2. (a) Geometrical and material properties of the cylindrical panel: R=2540mm, L=254mm, θ=0.1rad, thickness h=6.35mm, E=3102.75 MPa, ν=0.3 (b) Equilibrium curve of the perfect cylindrical panel.

Figure 3 .Figure 4 .Figure 5 .

 345 Figure 3. Projections and 3D representation of the fold line.

Figure 6 .

 6 Figure 6. Series improvement with Padé approximants.

Figure 7 .

 7 Figure 7. First buckling mode of the panel without imperfection.

Figure 8 .

 8 Figure 8. Projections and 3D representation of the "load-displacement-imperfection" fold line.

Figure 9 .cFigure 10 .

 910 Figure9. Fold line and equilibrium curves with associated limit points for different values of the imperfection amplitude (η = 0, η = 2.36h and η = -3.49h).

  nl (φ 0 ) + B T nl (φ 0 ) D B0 + G T Ψ0 G dV [25] F 1 = Ω BT 0 D B nl (u 0 ) + G T Ŝ0 G u ⋆ 0 dV [26] F 2 = Ω BT 0 D B nl (φ 0 ) + B T nl (φ 0 ) D B nl (u 0 ) + G T Ψ0 G u ⋆ 0 dV [27]Appendix I.2 : Thickness imperfection

  and Keller [KEE 73]. It has subsequently been studied or used by Moore and Spence [MOO 80], Weinitschke [WEI 85], Wriggers et al. [WRI 88] and Wriggers and Simo [WRI 90], among others, for the precise calculation of limit and simple bifurcation points in the case of perfect structures. Indeed, in this particular case (Λ = 0), the solution of [3] is restricted to a single point. The mathematical study of [3] with Λ varying has been carried out by Jepson and Spence [JEP 85] and it has been numerically investigated by Eriksson [ERI 94] and Deml and Wunderlich [DEM 97] for sensitivity analyses, in a finite element framework and using incremental-iterative strategies.

Produit par HAL -16 Sep 2014

Appendix I : Directional derivatives of F In this appendix, the procedure for obtaining the directional derivatives K T , K φ , F 1 and F 2 is presented for both shape and thickness imperfections.

In the case of geometrical nonlinear elasticity, [1] is cubic with respect to u. This cubic expression is not very suitable for asymptotic expansions. A quadratic expression is preferred. It is achieved by introducing the stress-strain relation as an additional equation.

[1] can thus be replaced by the following equivalent system

We have used the classical B operator defined by B(u)=B l + B nl (u) where B l and B nl (u) are the classical operators expressing the linear and nonlinear parts of the Green-Lagrange strain [CRI 97]. S is the second Piola Kirchhoff stress operator and D is the classical elasticity operator function. The first equation stands for equilibrium. It is quadratic with respect to the set of variables (u, S). The constitutive law has been introduced in order to make both equations quadratic. The additional equation [4] is then obtained by differentiation with respect to u and reads

where Ψ is the stress associated with the eigenvector φ.

Appendix I.1 : Geometrical shape imperfection

In order to get a scalar extra parameter, the imperfection is written as

where u ⋆ 0 is a fixed displacement which gives the shape of the imperfection and η is its amplitude. With these notations, the introduction of the imperfection transforms equations [15] and [16] into

Produit par HAL -16 Sep 2014 order to reduce the order of the equations and to turn them into a quadratic form. This is achieved as follows

[28]

Through the introduction of the variables S, γ, Ψ, ε and e, each equation of the previous systems is quadratic. By this way, the number of products which follow from the introduction of the series will be significantly reduced. Each of the variables u, S, γ, φ, Ψ, ε, D, e, h and λ is then expanded into a power series and replaced in [28] and [29], yielding at order 1

The successive elimination of all the additional variables in [30] is performed as follows. e 1 is first substituted into D 1 , then D 1 and γ 1 are replaced into S 1 and finally S 1 is reported into the equilibrium equation. Applying the same process to [31] yields again equations [22] and [23], but this time F 1 and F 2 read

where the notation D = (D 1 + 3h 2 0 D 3 ) has been introduced. The expressions of K T and K φ are the same as in appendix I.2 excepted that D must be replaced by D and B0 by B 0 because of the new initial conditions.