Three-level inventory deployment for a luxury watch company facing various perturbations

Abstract A well-known Swiss watch brand, active in the top-end luxury market, is facing a complex inventory deployment problem where watches of different models (more than 100 different models) must be dispatched first to wholesalers to finally reach the shops where consumers come in. Along the way, different perturbations are expected at three levels (production plan, demand, and dispatching), and accurate reactions must be taken to fit to these uncertainties. Solution methods are proposed to solve realistic instances. Results show the relevance of the methods and the robustness of the solutions.


Introduction
In this paper, a well-known Swiss watch brand (denoted SWB), which is active in the luxury industry, is facing a multi-level inventory deployment problem called (P), where watches (called items) of different models (called products) are produced, and must be sent to the final shops while satisfying various constraints. Inventory deployment is often regarded as a very complex problem by companies, as different constraints and perturbations often raise along the different levels (i.e. from the production to the final consumers), and the dispatching system must react accordingly. The supply chain is composed of three different levels, namely a single factory, the wholesalers (WS), and the final destinations, which are the shops where the goods are sold to the consumers. The factory is the only place where the goods are produced. For each considered geographical region, there exists at least one wholesaler owned by SWB which then delivers to its associated set of shops. Each shop can only be delivery by its single associated wholesaler. As far as SWB is concerned, the shops are the final points of interest, and the final consumer, who actually buys an item, is not integrated in the system as the shops are usually not owned by SWB. Thus, in the following, the customers of SWB are the WS, and the consumers are the persons actually buying the items at the shops' level. Unfortunately, in (P), three different perturbations are encountered along the supply chain: first at the production plan level, then at the dispatching level (WS), and finally on the demand (i.e. the actual demand is different from the forecasted demand). These perturbations are, respectively, due to the suppliers, the behaviours of the WS (on which SWB has not a total control), and the behaviours of the consumers.
The considered market involves a relatively small number of watches but a very high price per watch. Even if we tackle a three-echelon inventory deployment problem which would involve three types of decision-maker (i.e. plant, WS, shops), the only true and crucial decision-maker is located at the plant level, which has to decide on the deployment of the produced items (and thus also on the assortment planning), following an imposed push policy (from the plant to the shops, transiting through the WS). In other words, no inventory is kept at the production and WS levels (for safety, strategic, and confidential reasons). As soon as an item is produced, it has to be shipped to a shop (through a wholesaler), and the planned shop is decided at the production level. We, however, account for the behaviour's of the decision-makers located at the WS' level. Indeed, we propose to model the perturbations on their dispatching decisions (which do not always match the decisions made at the plant level). Therefore, the shops (resp. WS) do not order the watches to the WS (resp. plant), but they only provide them with information on their on-hand inventories. In other words, the WS and the shops have to agree on getting what the plant level has planned to send to them. Therefore, given this restrictive and imposed push approach (indeed, SWB only wants to improve its current practices), no inventory management policy is used along the supply chain (e.g. reorder point s with a fixed order quantity Q, reorder point s with order-up-to-level S).
As the watches produced by SWB are expensive, the suppliers have a key role because the components used to produce the luxury items are expensive and rare. Thus, there is no long list of suppliers which can provide the components. The number of products and the number of items yearly produced are both increasing with the demand, which makes the plant reach its production capacity limit. From a SWB perspective, the consumer wants his/her watch when he/she has decided to buy it. There are thus inventories in the stores, from which the final consumer is immediately delivered as soon as he/she has decided to buy a watch. As a consequence, a perturbation of the production plan due to suppliers has a strong impact on the availability of the items at the shops' level and is hard to manage, as every item that is not produced due to unreliable suppliers must be produced later and therefore disturbs the production plan of the following weeks.
The planning horizon is a year, and the time unit is a week. (P) is a three-level inventory deployment problem, where the main decision to make is the number of items of each product to send (from the plant and through the WS) to each shop on a per week basis. A solution is an aggregation of these numbers for a whole year. Even if the final clients of the plant are the WS, SWB must integrate the shops in their supply chain, as they are the only actors dealing with the consumers' demand. The objective function involves three different components: shortages, rarity, and inventory. The rarity of a solution is defined as the opposite of the diversity. For a given time, the diversity indicator of a shop is the number of different available products it contains. Then, the diversity of a solution is the summation of the diversity indicators over all the shops and time periods. In an ideal solution, at least one item of each product is represented in each shop for each time period. Thus, for a given time period, if a product is not represented in the assortment of a shop, then a rarity penalty is applied. One can already observe that rarity and inventory are conflicting components of the objective function. Indeed, a low rarity value leads to high inventory costs.
Most of the literature involves only two levels. There is no literature dealing with such a complex problem (P) with unreliable behaviours of suppliers, WS, and consumers. Moreover, the perturbation on the production plan can result in an infeasible solution, and thus, a specific repair procedure is requested. This paper aims at proposing dedicated solution methods and a simulator to improve the quality of the inventory deployment of SWB. The simulator is used to evaluate the quality of a solution according to the above described random events and is helpful to conduct a sensitivity analysis, allowing to identify the key perturbation parameters (i.e. the ones which have a significant impact on the objective function). One of the proposed solution methods shows that the use of relevant distance functions can be very powerful for tackling (P), where exact methods are not appropriate on realistic instances. Considering the computation time limit imposed by SWB, the best method is a matheuristic, which combines the accuracy of an exact method with the speed of a heuristic. If larger time limits are allowed, the use of local search approaches is beneficial.
The paper is organized as follows. In Section 2, an exhaustive literature review is provided. In Section 3, the three different perturbations faced by SWB are first exposed, and then, a formal description of (P) is proposed. Section 4 describes the different solution methods used to provide a solution S, and the simulator is designed to evaluate the quality of a solution S. The experiments are presented in Section 5. Section 6 concludes the paper and identifies possible future works.

Literature review
As explained in Glasmeier (1991), the Swiss industry lost control on the watch industry in the 1990s, due to the different actors gaining power in the technological field. But, as of today, Swiss luxury watches are selling in an increasing fashion each year (FSWIFH, 2014), which shows how the Swiss luxury industry was able to overcome its past drawbacks to return to expansion. Switzerland is not the only country facing luxury expansion. In Brun et al (2008), the authors state that in Italy, the luxury segment was worth 26 billion Euros per year in 2006 and is growing. The paper aims at investigating 12 luxury Italian firms of different luxury segments, to define the behaviours of each segment in the supply chain strategies. Extending the previous paper, the authors of Caniato et al (2011) propose to derive different clusters of companies and identify which supply chain strategy is applied in each cluster. In Catry (2003), a study is proposed on how luxury companies are adept of pretending to sell rare and exclusive products, even if the quantity of each product is so important that it does not belong to rarity anymore.
Usually, the actors of the Swiss watch industry started with a small family-sized company and gained world-class profits, focusing on the products and having only basic knowledge of supply chain management. This trend tends now to be reversed and many companies are trying to optimize their supply chain systems. In Hnaien et al (2010), a two-level assembly problem is tackled using a multi-objective approach based on genetic algorithms. Uncertainties appear in lead-times, and the main objective is to maximize the service level from the client point of view. In Silver and Zufferey (2005), a known demand and uncertain lead-times are expected, where the perturbations of the lead-times depend on the season. In such a situation, heuristics are proposed and their performance are analysed. Extending the previous paper, perturbations on the lead-times are again expected in Silver and Zufferey (2011), but the company faces a planned shutdown period. Simulation is used to assess the quality of the proposed solution methods.
As also discussed in the above papers, to correctly assess the quality of the models and methods, simulation has to be used. With regard to simulation, different interesting surveys can be found in Banks (1998), Law (1991), Terzi and Cavalieri (2004). A reverse logistics problem is tackled in Kara et al (2007), where defective products must be brought back to the factory from the client level. The simulator is used to assess the robustness of the presented model. A simulator is used in Daniel and Rajendran (2005) to tackle an inventory deployment problem using a genetic algorithm. The robustness of the algorithm is tested with different supply chain settings. In Miranda and Garrido (2004), an inventory control problem is considered with stochastic demand, using either a nonlinear mixed integer model and a heuristic based on Lagrangian relaxation. The exact model achieves impressive results on this problem, whereas the heuristic is used for larger instances only. In Ben- Daya and Hariga (2004), a problem with stochastic demand and dynamic lead-times is investigated and an exhaustive sensitivity analysis is performed for the important parameters (e.g. lot sizes, reorder points). In Schmitt et al (2010), a model is proposed for a problem with complete supplier disruptions. Different levels of uncertainty are modelled (e.g. demand, disruptions) and a closed-form approximation method is proposed. In Anupindi and Bassok (1999), the problem of supply contracts under quantity commitment and stochastic demand is studied. Simulationbased optimization approaches are proposed in Otamendi and Doncel (2012) for the entire transmission problem (i.e. arrival of ships, regasification, transportation, and distribution) faced by gas utilities companies and in Tsai and Zheng (2013) for a two-echelon inventory system with service-level constraints.
In Madadi et al (2010), a multi-level inventory management problem is tackled, where suppliers, warehouses, and retailers are considered, and transportation costs are encountered. Centralized and decentralized ordering models are compared. A two-level inventory problem with one warehouse and many retailers is proposed in Forsberg (1995), where the retailers face different compound Poisson demand processes, and the facilities apply order-up-to-S replenishment policies. A very similar problem is proposed in Axsäter (1993), where simple recursive procedures are proposed to evaluate the shortage costs. In Axsäter and Juntti (1996), distribution systems with stochastic demands are analysed using simulation. Advantages of echelon stock (Clark and Scarf, 1960) and installation stock (Lee and Moinzadeh, 1987;Moinzadeh and Lee, 1986) are assessed. In Segerstedt (1996), a mathematical formulation is presented for a capacity-constrained multi-stage inventory and production control problem. Numerical experiments for a trivial version of the problem are presented, and links among other problems are proposed [e.g. Kanban (Bard and Golany, 1991), cover-time planning (Segerstedt, 1991)]. In Yao et al (2009), three different inventory strategies are designed for a two-level inventory deployment problem, and managerial insights are provided. In Svoronos and Zipkin (1988), the authors propose methods to assess the performance of different multi-level inventory strategies. In Kunnumkala and Topaloglu (2011), the authors consider an inventory-distribution system consisting of one warehouse and multiple retailers. The retailers face random demand and are supplied by the warehouse. The warehouse replenishes its stock from an external supplier, and the objective is to minimize the total expected replenishment, holding and backlogging cost over a finite planning horizon. Linear-programming-based decomposition methods are proposed to tackle this problem. In Schildbach and Morari (2016), a scenario-based model is designed for multi-echelon supply chain management. The authors showed that except for a few special cases, optimal solutions are computationally intractable for systems of realistic size. Their approach can handle supply chains with stochastic planning uncertainties from various sources (e.g. demands, lead-times, prices).
In Melo et al (2009), an interesting survey is proposed on facility location and supply chain management, in addition to recent advances in optimization techniques for these problems and for inventory deployment. In line with the current trend but in contrast with the practice of SWB, an integrated inventory management model is depicted in Shafieezadeh and Sadegheih (2014) for multi-item multi-echelon supply chains. The benefit of managing jointly different decision levels will be highlighted in the sensitivity analysis performed in Section 5.3.

Presentation of problem (P)
In this section, the different perturbations are explained, followed by a formal problem description of (P).

Perturbations
As explained earlier, the complex process of producing watches suffers multiple uncertainties. Three different perturbations are exposed, namely perturbations on the demand, perturbations on the production plan, and finally perturbations on the dispatching of the items at the WS' level. The intervals in which the perturbations parameters belong to are set thanks to the SWB knowledge of their supply chain.
3.1.1. Perturbations on the demand There is obviously a gap between the forecasted and the actual demands encountered at the shops' level. For SWB, the consumers are wealthy persons, who wants the item at the time they decide to buy it. From forecasted demands (given by SWB) and following the historical observations of SWB, we generate perturbations based on a normal distribution. Thus, each shop has its own mean and standard deviation for each product, used to generate the perturbation on the demand for each week. The means and standard deviations are set with the help of SWB to accurately reflect its reality.
3.1.2. Perturbations on the production plan For various reasons (suppliers being unreliable, problems in the deliveries of the components, etc.), two different types of perturbation on the production plan can occur every week. The first one sees some products being not produced at all (because the associated components were not delivered), whereas the second sees some items of some products being not produced (as the production was slower than expected). In other words, the production plan is decided (i.e. fixed), but it can suffer from perturbations. The challenge is thus to develop an inventory deployment method which can react appropriately to the uncertainty associated with the production plan.
Letp i t be the planned number of items of product i to produce at week t, and p i t be the corresponding actual number (i.e. with the perturbation). Each week, a maximum of d 1 per cent of models are not produced, due to non-deliveries from suppliers (i.e. p i t ¼ 0 for these models, instead of p i t ¼p i t ). The perturbation affects 2 weeks, and the corresponding production must be performed later. If a component is not delivered to the plant at period t, SWB assumes that it will be necessarily delivered by period t þ 2. During period t þ 2, SWB is then always able to adjust its production capacity in order to recover from this previous components' shortage. Thus, p i t ¼ p i tþ1 ¼ 0, and from a practical standpoint, we can always set p i tþ2 ¼p i tþ2 þp i t þp i tþ1 . SWB proposes that d 1 belongs to interval [5,15].
On the remaining 1 À d 1 per cent of the models (i.e. the ones which are produced), each model has d 2 per cent of risk to be slowed down. Slowed down means that at most d 3 per cent ofp i t is not produced at week t but at week t þ 1. At the end of week t, the total number of items not produced due to slowed down performances is Q t ¼P t À P t , whereP t ¼ P ip i t and P t ¼ P i p i t . On the total numberP t of items expected to be produced at week t, a given threshold of e per cent of items are allowed to be slowed down, meaning that if Q t e ÁP t , then the factory does not take any action. But if Q t [ e ÁP t , then C t ¼ Q t À e ÁP t watches must be compensated (meaning that items of different models are produced instead of planned ones) with the following rule: models which are not slowed down at week t and planned at week t þ 1 are eligible for compensation, and the models which were not expected to be produced at week t (i.e.p i t ¼ 0), but are expected to be produced at week t þ 1, can be eligible for compensation as well. For the eligible models, some watches (randomly selected, one at a time) have to be produced at week t, until quantity C t is fulfilled, as explained further in Section 4.2. SWB proposes the following intervals for the above discussed parameters: d 2 2 ½45; 55, d 3 2 ½5; 15, e 2 ½5; 15. Note that, for each product, SWB has a fixed annual quantity to produce, which has to be met in any case. The current practices of SWB allow to assume that SWB is always able to adjust its production capacity if a compensation action is triggered.
3.1.3. Perturbations on the dispatching from the wholesalers As shown in Figure 1, the supply chain network linking the factory, the different WS, and the shops is a threelevel inventory model. The WS belong to SWB and are the only customers regarding SWB. It is expected that the different WS manage the received inventory and deliver items to their corresponding shops in a reliable manner, but unfortunately perturbations occur. Indeed, each wholesaler has its ''confidential'' list of preferred shops, and it can decide to ship an item to a preferred shop instead of the one initially planned by SWB at the plant level. The perturbation which can happen is the following: SWB estimates a target inventory for each final shop and sends the corresponding items to the associated wholesaler. Each shop j has a priority w j and let w max be the largest possible priority (i.e. the priority of the most important shops). Any wholesaler can perturb the SWB estimations by sending the items dedicated to a given shop j to a different shop j 0 . SWB assumes that d 4 per cent of the items does not go to the expected shop. On these d 4 , d 5 per cent would go to a shop with a higher priority. SWB proposes to use d 4 2 ½30; 50 and d 5 2 ½60; 80.

Formal description
To deeply understand (P), we propose a formal description. N models (i.e. products) are produced in the factory. The number of weeks in a year is W. Each week t, p i t items of model i arrive to the different WS around the world. At the end of each week, because of safety/strategic/confidential reasons, all the produced items are sent to the shops (through their associated WS), and thus, none are kept in inventory at the plant and WS levels. Such a push policy is imposed by SWB which wants to improve its current practices. As a consequence, no inventory management policy is used by the shops and the WS. J shops must be supplied by the different WS, and each shop j is associated with exactly one wholesaler m j (thus, we can have be the on-hand inventory of wholesaler m j regarding model i at week t. Each wholesaler delivers to its shops and tries to empty its inventory. Indeed, the WS are owned by SWB, which naturally prefers to put the inventory costs at the shops' level. Letx i;j t be a decision variable (managed by the only key decision-maker, located at the plant level) describing the number of items of model i supposed to arrive at shop j at week t from the corresponding wholesaler m j (it thus corresponds to a planned shipment from the plant). One can observe that thex i;j t 's and thep i t 's are related. Indeed, the former variable is a part of the latter, but two different times are involved as the shipment from the plant to the shop requires a lead-time. As explained in Section 3.1.3, each variablex i;j t suffers a perturbation, and let x i;j t be the resulting variable involving the perturbation. In other words,x i;j t (resp. x i;j t ) is an expected (resp. actual or observed) value. Moreover, a lead-time of L m j weeks occurs from the plant to the wholesaler m j . Another lead-time of L j weeks is then required to reach the shop j from the corresponding wholesaler m j . For each week t, a demandd i;j t is forecasted for model i at shop j and is subject to the perturbation depicted in Section 3.1.1. Thus, let d i;j t be the corresponding actual demand. As soon as the items x i;j t reach the shop j, the corresponding on-hand inventory of shop j is updated: . Thus, the inventory at period t is set as the inventory at period t À 1 augmented by the number of delivered items and decreased by the demand at period The overall objective function f to minimize involves three different components f 1 , f 2 and f 3 to minimize in a lexicographic order (i.e. no deterioration on f i can be compensated by improvements on f iþ1 ), as motivated in SWB (2014). More precisely, one can define f as c are chosen to respect the lexicographic order (which is ensured for example with a = 1,000,000,000, b = 1,000,000,000 and c ¼ 1). In preliminary experiments, we have encountered several solutions with equivalent values for f i while satisfying an upper bound on f iÀ1 (for i 2 f2; 3g). This actually means that all the components of f are relevant and none can be ignored. In other words, even f 3 can be the discriminating component. This observation is in line with the findings of Respen et al (2016), where lexicographic optimization (involving three components as well) was applied in an industrial context arising in the car industry.
Considering (P), f 1 is the expected shortage penalty for each item i for each week t in each shop j, which occurs if The last objective f 3 is described in Equation (3), which is a weighted inventory penalty.
Observe that the inventory penalty ðw max þ 1Þ À w j depends on the priority of the involved shop. Indeed, as SWB makes more margins with the shops with higher priorities, SWB prefers to allow more stock in such important shops. For this reason, the inventory penalty is lower for the shops with a higher priority.

Solution methods and simulator
To evaluate the actual value of a solution, it is mandatory to develop a simulator. Its input data are a solution, the production plan, the initial inventories (for each shop and each wholesaler, for each product), and the forecasted demand. The overall approach is presented in Figure 2. This decision process is in line with the decision-making literature, for which the reader is referred to Malakooti (2012) for having a general reference. A solver is first used to generate a solution S, which is then evaluated with the simulator. Such information (i.e. the solution S and its simulation along with the various objective function components) is then provided to the decision-maker, which can then decide to manually modify the solution (based on his/her expertise and the non-modelled elements) and simulated it again. The solver is made of a solution method, which can involve or not simulation. In the former case, the solver accounts for uncertainty (it will be for example the case for the proposed methods denoted BSS and BSSC). Two computation time limits are considered in this study: T ¼ 30 min and T ¼ 10 h. A strict time limit of T ¼ 30 min (including simulation) is imposed by SWB. This forbids the use of time-consuming methods such as metaheuristics. This restrictive time limit is imposed because of managerial reasons. For instance, after the use of the solution method, the involved decision-maker has to take the provided solution as an input and rework on it (relying on non-modelled features) based on the process depicted in Figure 2. Such an optimization process involving the decision-maker was also confirmed in other studies [e.g. Syntetos et al (2011)]. We will see in Section 4.3 that simulating one solution requires roughly 3 min on the used computer. For these reasons, constructive methods are appropriate to tackle (P) if T ¼ 30 min. However, in order to further evaluate the quality of the proposed constructive algorithms, a time limit of T ¼ 10 h is also considered, allowing the use of more advanced metaheuristics, namely a descent local search DLS and a tabu search TS.
In the remaining part of this section, three constructive solution methods are proposed in Section 4.1, the simulator is presented in Section 4.2, important computing time aspects are discussed in Section 4.3, and metaheuristics are proposed in Section 4.4 (for the case T ¼ 10 h).

Constructive solution methods
We first design three solution methods, which rely on a population of pop (parameter) solutions and use a constructive heuristic called SmartPriorityGreedy (denoted SPG) as an underlying method. SPG is depicted in Algorithm 1, where A i;w t is a set that contains all the shops with priority w that have a strictly positive forecasted demand for week t and model i, and A contains all the shops with the highest priority w max (regardless of the demand). Note that for each shop and each wholesaler, the initial associated inventory level (i.e. the available quantity) of each model is a given as input. Next, for each week t and each model i, the available items are the sum of the produced items in the factory and the items available at the WS' level. At each step, the algorithm attributes one item to a shop (through its associated single wholesaler). After intensive tests, the most efficient method (according to speed and results' quality) consists in randomly selecting the shop to which an item is assigned in the considered set A i;w t . Using more advanced methods (using f 1 , f 2 and f 3 for example) dramatically slows down the overall process, without improving the quality of the solutions. As an example, on an instance with N ¼ 143, W ¼ 47, J ¼ 192 and M ¼ 9, a method which computes f 1 , f 2 and f 3 at each step to attribute the best item to the best shop requires more than 15 h to return a solution! Here, the efficiency of SPG relies on the design and efficient use of the A i;w t 's. f 1 and f 3 are very influenced by the weights of the shops, and creating the A i;w t 's allows to very quickly and efficiently improve these objective components. If all the expected demands of model i of all the shops have been fulfilled for the considered week t, the last phase of SPG is triggered [see the last part of step (3c) in Algorithm 1]. This phase consists in shipping the remaining items (i.e. for which there is no demand) to the shops with priority w max : The first solution method, called BestSolution (denoted BS), generates a set of pop solutions using SPG and returns the best one (according to f). Following this method, BestSolutionSimulated (denoted BSS) generates a set of pop solutions with SPG, simulates the k (parameter) best ones (according to f), and finally returns the solution with the best simulated value. Finally, BestSolutionSimulatedClique (denoted BSSC) is a matheuristic (Boschetti et al, 2009).
To properly understand BSSC, let us first introduce a distance function distðŜ 1 ;Ŝ 2 Þ which returns the distance between two solutionsŜ 1 andŜ 2 . We propose to set distðŜ 1 ;Ŝ 2 Þ ¼ P t;i;j w j Á k i;j t , assuming thatx i;j t ðŜÞ is the value ofx i;j t in solutionŜ, and k i;j t is set as in Equation (4). Therefore, this distance function assumes that two solution components are distant if one is zero and the other is strictly positive (regardless of the amplitude of the positiveness), which is relevant in the context of SWB, which wants at least one item of each product available in each shop.
Let G ¼ ðV; EÞ be a graph with vertex set V and edge set E. A clique C of G is a subset of V such that its vertices are pairwise adjacent. In addition, let wðv; v 0 Þ be the weight associated with edge ðv; v 0 Þ 2 E. We define here the weight w(C) of a clique C by 1 2 Á P v;v 0 2C wðv; v 0 Þ. The coefficient 1 2 ensures that the weights are not counted twice. Further, the kclique with maximum weight problem consists in searching the clique C of size k in G which maximizes w(C). This problem is also known in the literature as the maximum diversity problem, for which extensive computational experiments are performed in Marti et al (2013) to compare 10 heuristics and 20 metaheuristics. Surveys on graph theory and cliques can be found in Brandstädt and Spinrad (1999), Pardalos and Xue (1994), Schaeffer (2007). This problem can be exactly solved with an integer linear program (ILP) using CPLEX. Let The ILP model is given in Equations (5-10). Constraints (6) and (7) are used to compute correctly the value of y v;v 0 , and Constraints (8) ensure that the clique size is k. x X x v 2 f0; 1g 8v ð9Þ BSSC starts by generating a set of pop solutions with SPG, then builds a set B with the best b (parameter) generated solutions, and on this set computes distance between each pair of solutions. The above ILP clique problem is then solved to optimality with CPLEX [a vertex v represents a solution S v 2 B and wðv; v 0 Þ ¼ distðS v ; S v 0 Þ]. BSSC finally simulates only the k solutions associated with the provided clique. At the end, BSSC returns the solution with the best simulated value. The clique approach used in BSSC allows to select the solutions that differ the most (among the most promising ones), with respect to their structures. In addition to being efficient, this technique allows the decision-makers to select a solution among highly different solutions.

Simulator
The simulator returns a value f(S) of a solution S after having added the different perturbations to the corresponding solution without perturbationŜ (as explained in Section 3.1) in that order: generate the demand perturbations, perturb the production plan, and finally dispatch. Algorithm 2 proposes a sketch of the different simulation steps. After the second step, which perturbs the production plan, the solutionŜ could possibly not be feasible anymore, as thep i t 's could be different from the p i t 's (i.e. some items are unfortunately not produced, and others are produced in counterpart). Thus, a repairing step is performed, which ensures that the resulting solution remains feasible. Let f j k be the contribution of shop j to the objective function component f k . To repair a solution, is computed for each shop j, and the resulting f j 's are sorted in an increasing fashion. If for a given week, ax i;j t has to be reduced (because of the production plan perturbations, some watches do not exist anymore), then they are one by one removed from the shops with the smallest f j value. On the opposite, if some items must be added (because of the compensation production), then they are added one by one to the shops with the largest f j value. A pseudo-code of the repair procedure is given in Algorithm 3. Note that even after the perturbations on the dispatching, solution S remains feasible because the watches are managed.
At the very end, the simulator returns a value f(S) for the input solutionŜ. To properly capture the operating mode of the simulator, let us propose a small example. Assume that an instance involves one product, 2 weeks, one wholesaler, and two shops (with w 1 ¼ 2, w 2 ¼ 1, and thus, w max ¼ 2). All initial inventories are empty, and all the lead-times are zero. The production plan isp 1 ¼ ð2; 3Þ and the forecasted demand d 1;1 ¼ ð1; 3Þ andd 1;2 ¼ ð1; 1Þ. The solver, based on the expected production plan and demand, proposes a solutionŜ withx 1;1 ¼ ð1; 3Þ andx 1;2 ¼ ð1; 0Þ, as it fulfils in priority the demand for the first shop. Then, the actual demand is generated as d 1;1 ¼ ð1; 2Þ and d 1;2 ¼ ð2; 1Þ, whereas the actual production plan is generated as p 1 ¼ ð3; 3Þ (we produce more than expected, due to a perturbation which occurred before the current events). The solution then needs to be repaired, and the system returns a modified solutionx 1;1 ¼ ð1; 2Þ and x 1;2 ¼ ð2; 1Þ, which is optimal. But the behaviour of the wholesaler perturbs the solution such that x 1;1 ¼ ð1; 3Þ and x 1;2 ¼ ð2; 0Þ. Finally the objective function is computed as f ðSÞ ¼ 1 þ 1 4 þ 1 ¼ 9 4 (using a ¼ b ¼ c ¼ 1). In order to be robust, sim (parameter) runs of the above described simulator are performed with solutionŜ as input. The value f(S) is then the average value over the sim runs. Let f sim ðSÞ be the average value of f(S) after sim runs of the simulator. Parameter sim was tuned such that f sim ðSÞ and f simÀ1 ðSÞ differ by \1. After intensive tests on various instances, we have set sim ¼ 40. More generally, if sim \40, there exists minor variations of the average value of f(S) over the sim runs, whereas if sim ! 40, only very negligible variations of f sim ðSÞ are observed.

Computing time considerations
Remind that the time limit T allowed by SWB is 30 min. A typical instance of SWB contains 47 weeks, from 7 to 9 WS, from 160 to 250 shops, and from 100 to 200 models (these intervals are voluntarily consequent, to respect a SWB nondisclosure agreement). SPG requires half a second to get a solution. A complete simulation (i.e. sim = 40 runs of the simulator on one solution) requires around 3 min. Based on such constraining time considerations, preliminary experiments lead to the following parameter setting: pop ¼ 100, b ¼ 25, and k ¼ 5. With such values, CPLEX requires about 5 min to optimally solve the ILP clique problem [see Equations (5)-(10)]. BSS and BSSC need between 20 and 30 min to complete, whereas BS requires only 15 s. If one would reduce the computing time of BSSC, an option is to solve the ILP clique problem with a (meta)heuristic [as proposed in Marti et al (2013), but without any guarantee on optimality] instead of with CPLEX. However, as BSSC meets the 30-min time limit imposed by SWB, we propose to keep the use of CPLEX.

Local search methods
Assume that an objective function f has to be minimized within a given time limit T. A local search algorithm starts from an initial solution (usually provided by a constructive heuristic) and tries to improve it iteratively. In each iteration, a move (i.e. a slight modification) is performed in order to transform the current solution into a neighbour solution. In a descent local search (DLS), the best move (according to the reduction of f) is performed in each iteration, and the process is restarted from another initial solution as soon as a local minimum is found (the best encountered solution is provided at the end). In contrast, tabu search (TS) continues performing moves until the time limit T is reached, but in order to prevent cycling (i.e. coming back to an already visited solution, like the last encountered local minimum), the reverse of the most recent moves are forbidden (they are called tabu moves) for tab (parameter) iterations. It results that in each iteration of a basic TS, the best non-tabu move is performed. The reader is referred to Gendreau and Potvin (2010) for a book on metaheuristics, whereas Zufferey (2012a) proposes general guidelines on metaheuristics to efficiently design them depending on the encountered problem.
Considering (P), remember that each decision variablex i;j t is the number of items of model i supposed to arrive at shop j at week t from the corresponding wholesaler m j (it corresponds to a planned shipment from the plant). For fixed values of t and i, a straightforward move (both for DLS and TS) consists in shipping an item of model i to another shop j 0 instead of the currently planned shop j, while possibly involving another week t 0 6 ¼ t because of lead-time considerations. In other words,x i;j t (resp.x i;j 0 t 0 ) is reduced (resp. augmented) by one unit. As discussed in Section 3.2, as long as a model i is not represented in a shop j, a shortage penalty is applied by SWB. In order to focus on the most promising moves, the following two conditions have to be satisfied:x i;j t ! 2 andx i;j 0 t 0 1. In other words, for model i, an interesting move should not lead to a shortage penalty in shop j, it should not augment too much the inventory level in shop j 0 , it should have a chance to improve the assortment planning of shop j 0 , and it should have a chance to remove a shortage in shop j 0 . When a move is performed in TS, it is then forbidden to augment (resp. reduce) againx i;j t (resp.x i;j 0 t 0 ) for tab iterations, where tab is an integer randomly generated in interval [5,15] (thus, each tabu feature has its own tabu duration). Preliminary experiments have confirmed that such a way of managing the tabu status is more efficient than having a single value of tab for the whole search process.
We have two evaluation measures: (1) the expected cost f ðŜÞ of the expected solutionŜ; (2) the simulated cost f ðsimÞ ðSÞ of the actual solution S associated withŜ, where sim is the number of replications (tuned to 40 for having reliable and robust results, as discussed in Section 4.2). As the simulation duration is significant (3 min/solution), not all the neighbour solutions can be accurately simulated in each iteration (otherwise \5 iterations would be performed in 10 h). For each iteration, a sensitive issue is thus the way to build a small but promising set of neighbour solutions to evaluate with f ðsim¼40Þ . The following process is proposed: (1) generate a random set of N 1 (parameter) neighbour solutions of the current solutionŜ; (2) evaluate the best (according to f and using incremental computation) N 2 N 1 solutions with f ðsim¼1Þ ; (3) evaluate the best [according to f ðsim¼1Þ ] N 3 N 2 solutions with f ðsim¼40Þ ; (4) the selected neighbour solutionŜ 0 ofŜ is finally the best [according to f ðsim¼40Þ ] of these N 3 solutions. In other words, the best expected solutions are first partially (but quickly) simulated, and only the resulting most promising solutions are then accurately evaluated with a full simulation. After a preliminary set of experiments, parameters ðN 1 ; N 2 ; N 3 Þ were, respectively, tuned to (20, 10, 3). With such a parameter tuning, a single iteration requires roughly 10 min of computation.
In order to generate an initial solution, both DLS and TS use BSSC. However, anytime DLS is restarted, the new initial solution is provided by BS, otherwise BSSC (which requires about 30 min) will use too much of the time budget T allocated to DLS.

Results
Tests were performed on an Intel Quad-core i7 @ 3.4 GHz with 8 GB DDR3 of RAM memory. Section 5.1 describes the instances used to evaluate the quality of the different solution methods. Section 5.2 shows the experiments on the realistic instances (considering the imposed time limit of T ¼ 30 min). Section 5.3 conducts a sensitivity analysis of the perturbation parameters. Finally, Section 5.4 presents the results obtained with DLS and TS when the time limit is extended to T ¼ 10 h.

Generation of the instances
Based on the various data provided by SWB, we have developed an instance generator, which is able to quickly generate some realistic instances. An instance is composed of the following parameters: the number N of models, the number W of weeks in a year, the number J of shops, the number M of WS, the perturbation parameters ðd 1 , d 2 , d 3 , d 4 , d 5 and eÞ, the production plan, the forecasted demand, the lead-times (from the factory to the WS, and from the WS to the shops), the shop priorities, the links WS-shops, the inventory levels at the WS and at the shops. The data parameters were reasonably fixed after discussions with SWB. In the reference instance, we use We have generated a total of 64 instances around the reference instance by varying d 1 , d 2 , d 3 , d 4 , d 5 and e.

Results on realistic instances
We propose below a summary of the results on the realistic instances, whereas the detailed results (and their analysis) can be found in the Appendix.
Aggregated results are presented in Table 1, where for each triplet ðd 1 ; d 2 ; d 4 Þ, we have generated eight instances by considering different values for the parameters ðd 3 ; d 5 ; eÞ. More precisely, d 3 2 f5; 10g, d 5 2 f60; 70g and e 2 f5; 10g. For each instance, let f H be the best value ever found for a solution. For each method and each instance group (i.e. for each line of Table 1), average gaps with respect to the average value of f H are provided. We decided to focus on the parameters ðd 1 ; d 2 ; d 4 Þ because they are the roots of the perturbations. More precisely, d 1 is the key perturbation parameter for the non-production perturbation, d 2 for the slowed down perturbation, and d 4 for the dispatching perturbation. We can roughly observe that both BSS and BSSC outperform BS on most of the instance groups, whereas BSSC gets the best results on average, which highlights that the proposed combination of a constructive heuristic with an exact model (based on an appropriate distance function) seems to be relevant for (P). Interestingly, when BSSC obtains the smallest gap on f 1 , it usually reaches the lowest gap on f 3 . On the opposite, when BSS gets the smallest gap on f 1 , it generally reaches the lowest gap on f 2 but not on f 3 . This tends to prove that BSSC is better at optimizing both f 1 and f 3 . Note that all the gaps for f 2 are relatively small (as also confirmed in the Appendix). Regarding the number of times that BSS or BSSC reach the lowest gaps, both methods are comparable. But when BSSC gets the best results on f 1 and f 3 , the corresponding gaps are significant when compared to the other two methods. Table 2 summarizes (over all the 64 instances) the performance of the different methods on the three objectives. The second row (labelled with ''Number of zero gaps'') gives the number of times that the corresponding algorithm founds the best value (over all the tests with all the methods). The last row shows the average percentage gap of each method with respect to the average value of f H . These results show that BSSC is the most efficient method on f 1 and f 3 , but it can be outperformed on f 2 , which is not surprising because of the lexicographic approach involving conflicting functions.

Sensitivity analysis
In this subsection, a sensitivity analysis of the most sensitive perturbation parameters is conducted. Figure 3 shows a wider variation of parameters d 1 , d 4 and d 5 (respectively). For each part of the figure, the other involved parameters are set to the reference values. The sensitivity analysis of the other parameters (i.e. d 2 ; d 3 ; e) was not very conclusive and therefore not presented here. Figure 3 (left part) shows that d 1 has an impact on both f 1 and f 2 (as f 1 decreases by almost 1% if d 1 decreases by 4%). Finding more reliable suppliers could lead to such a d 1 reduction. Figure 3 (middle and right parts) shows that both d 4 and d 5 are the most sensitive parameters, as they significantly reduce f 1 and f 3 (even if in the middle part of the figure, it leads to a small augmentation of f 2 ). This is due to the WS, which are likely to sacrifice the considered shop and favour higher priority shops. More precisely, as the shops with the lowest priority are the most present, the watches will rather reach shops of priority 2 or 3 and therefore decrease f 1 and f 3 (as shortages are less penalized if more watches arrive, and shops of higher priorities are less penalized on inventories). Therefore, f 2 increases for the shops of lower priorities, which are more present. These conclusions, regarding the WS, were acknowledged by SWB. This analysis is in line with the recent literature and trends, where in Yu et al (2012), a consignment policy is proposed for inventory management in supply chains, based on strong interaction and reliable collaboration between the vendor and buyer levels.
To measure the impact of the perturbation parameters on a solutionŜ, we propose to use DistðŜ; SÞ ¼ P t;i;j k i;j t , with k i;j t defined as in Section 4.1. Function Dist measures the structural difference between the expected solutionŜ and the actual solution S (after simulation). To better capture the impact of the perturbations, we propose to compute such a distance value for each priority level. More precisely, at the end of each single run of the simulation, the distance difference (in %) betweenŜ and the returned actual solution is computed for each priority level (i.e. each gap is divided by where J w is the number of shops with priority w). These gaps are then averaged over the sim runs. Such a test was performed 100 times on two different instances, and the results are averaged for each priority level. The used instances are ðd 1 ; d 2 ; d 3 ; d 4 ; d 5 ; eÞ 2 fð5; 40; 5; 30; 60; 5Þ, ð10; 50; 10; 40; 70; 10Þg (i.e. the instances with the smallest and largest perturbations, respectively). The resulting gaps are presented with the notation (w 1 ,w 2 ,w 3 ), where w i is the gap corresponding to the shops with priority i. The obtained gaps are (2.92, 3.11, 3.69%) for the least perturbed instance and (3.15, 3.50, 4.49%) for the reference instance. As an example, the shops with priority 3 have a gap of 4.49% on the reference instance. The results show that the larger gaps concern the most important shops, which is consistent as such shops are the least represented in the distribution network.

Results with T = 10 h
Three methods are compared in Table 3, namely BSSC, DLS, and TS. As shown in Table 1, aggregated but representative  On the other hand, as TS outperforms DLS (it is the case for each instance group), it seems more appropriate to employ tabu status rather than restarts in order to better control and diversify the search process.

Conclusion and future works
In this paper, we investigate a relevant and complex problem (P) faced by a well-known Swiss watch brand. We propose powerful solution methods to tackle (P) and showed the superiority of the proposed matheuristic compared with standard heuristics when considering the strict computation time limit imposed by the company. Moreover, we highlight the most sensitive parameters on which luxury companies should focus on. If larger computation time limits are allowed, it was showed that local search approaches (e.g. tabu search) can improve the results by above 2% on some objectives. Future works include developing other metaheuristics for the problem [e.g. among the different existing ant algorithms (Zufferey, 2012b)] or modelling a more complete problem that integrates more types of perturbations. One can also work on improving the speed of the proposed algorithms by jointly using lower bounds and solution space reduction techniques [as performed in Hertz et al (2005)]. Moreover, SWB could evaluate and optimize an integrated supply chain where suppliers, production, and dispatching would be centralized in a common system. This would allow SWB to reduce the perturbation parameters by gaining control over the complete supply chain. Appendix: Results on the realistic instances Table A1 shows the results for the three solution methods when d 1 , d 2 , d 3 , d 4 , d 5 and e vary. We have set r ¼ 2 (which is the standard deviation of the demand perturbation), as larger values are totally not realistic for SWB. Two different values are presented for each parameter, the first is lower than the reference value, and the second is the reference value. The goal is to compare the performances of the algorithms and the impact of lower perturbations on the objectives. The first six columns give the values of each parameter. The seventh column, called f H 1 , is the best f 1 value returned by one of the algorithms on the considered instance. The same applies for the next two columns regarding f 2 and f 3 . Starting at column ten, the next three columns show the gaps between the f H i 's and the f i 's for method BS. The same information is then given for BSS and BSSC. The last two rows show the number of times that the corresponding algorithm found a gap of zero, followed by the average values of the different gaps. Note that the first row is the instance with the lowest perturbations, and the last row shows the instance with the largest perturbations (which is the reference instance). The results shows that BSSC outperforms BSS on both f 1 and f 3 , whereas BS is significantly overpowered by both BSS and BSSC. Regarding f 2 , BSS finds zero gaps 29 times, versus 15 times for BSSC, but the gaps for BSSC are so small when BSS finds zero gaps that no conclusion can come out. Moreover, as f 1 has the priority over f 2 , it is not surprising that if BSSC outperforms BSS on f 1 , then BSS can sometimes outperform BSSC on f 2 . When comparing the last row, the averages over the complete set of instances clearly show that BSSC is the most efficient method, as it shows the lower percentages for f 1 and f 3 , and very low gaps for f 2 .
Regarding the impact of a single parameter variation (for example d 1 2 f5; 10g), and leaving the other parameters unchanged, we can see that the impacts on f 1 and f 3 are significant, but not on f 2 . The reasons to this are detailed in Section 5.3. The parameters that induce the greater reductions on f 1 and f 3 are obviously d 4 and d 5 . Both parameters are used when encountering the perturbations involved by the WS' behaviours. The improvements are such that it would be very important for SWB to reduce the perturbations by either finding more reliable suppliers (for d 1 ; d 2 ; d 3 ) or by having a better control on the WS (for d 4 ; d 5 ). To properly assess the impact of such reductions, Section 5.3 proposes a more detailed sensitivity analysis.