{'fraction_non_alphanumeric': 0.029875, 'fraction_numerical': 0.012964285714285715, 'mean_word_length': 4.309661515122784, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Mobile Production Vehicle Routing Problem: Using 3D Printing in Last Mile Distribution

We study a new variant of the vehicle routing problem, called the Mobile Production Vehicle Routing Problem (MoP-VRP). In this problem, vehicles are equipped with 3D printers, and production takes place on the way to the customer. The objective is to minimize the weighted cost incurred by travel and delay of service. We formulate a Mixed Integer Programming (MIP) model and develop an Adaptive Large Neighbourhood Search (ALNS) heuristic for this problem. To show the advantage of mobile production, we compare the problem with the Central Production Vehicle Routing Problem (CP-VRP), where production takes place in a central depot. We also propose an efficient ALNS for the CP-VRP. We generate benchmark instances based on Vehicle Routing Problem with Time Windows (VRPTW) benchmark instances, and realistic instances based on real-life data provided by the Danish Company 3D Printhuset. Overall, the proposed ALNS for both problems are efficient, and we solve instances up to 200 customers within a short computational time. We test different scenarios with varying numbers of machines in each vehicle, as well as different production time. The results show that these are the key factors that influence travel and delay costs. The key advantage of mobile production is flexibility: it can shorten the time span from the start of production to the delivery of products, and at the same time lower delivery costs. Moreover, long-term cost estimations show that this technology has low operation costs and thus is feasible in real life practice.


Introduction
In the modern supply chain system, balancing the interests of manufacturers, distributors, and customers has always been a difficult problem. Manufactures and Distributors want to minimize their operation costs.
Customers, meanwhile want to receive a good quality of service. Among the various costly factors behind the operations, inventory management is a long-time problem that has complicated the production and circulation of commodities. Annual reports at Amazon, for example, have listed the inventory problem as a risk factor within the operational process for two consecutive years (Amazon, 2017;Amazon, 2018). In 2017, Amazon reported that a 1 percent inventory increase would cost an additional 180 million USD (Amazon, 2017). With the rapid development of industrial manufacturing, an effective solution called make-to-order (MTO) was proposed to reduce such inventory costs.
In MTO mode, when the central warehouse receives an order, production begins there directly. The number of goods to be produced is the same as the number of orders received (i.e., no excess good is made).
When production is completed, trucks deliver the assigned goods to the customers, as shown in Figure 1. To minimize the total cost, the practitioners want to find the optimal production scheduling and transportation scheme. This need has raised an optimization problem, which we have named the Central Production-Vehicle Routing Problem (CP-VRP).
Although the on-demand central production in MTO mode can eliminate inventory costs, the potential drawback is that it is not flexible enough because waiting for the completion of production may delay delivery to the customers. This disadvantage does not conform to customers' wish for a high-quality service, which may ultimately lead to a decline in both service quality and brand reputation.
To improve and innovate the logistics model, Amazon applied for a patent (US Patent 9684919) in 2015 (Amazon Technologies, 2015;Boyle, 2018) that combines 3D printing technology and MTO mode. With this patent, a car-mounted 3D printer produces the products on the way to the customer, as shown in Figure   2. By applying this technique, production and transportation are synchronized and delivery times can be largely reduced. This zero-inventory, synchronous production, and transportation technology can meet people's demand for customization and prompt delivery, all while reducing the cost of logistics. Thus, applying this technology can be a good strategy to balance the interests of manufacturers, distributors, and customers in the modern supply chain system. With the further upgrading of industrial manufacturing and the continuous maturation of 3D printing technology, the era of intelligent manufacturing and personalized customization is coming.
Thus, this new technology will have great potential for large-scale use in the future.
The integration of the car-mounted 3D-printing factory and the vehicle routing problem is denoted as the Mobile Production Vehicle Routing Problem (MoP-VRP). Unlike the CP-VRP, in which delivery starts only after production is completed at the depot, production and distribution in the MoP-VRP take place synchronously. This setting makes the MoP-VRP more complex than the CP-VRP, as production and delivery are more closely connected. This close connection also makes the existing methods proposed for the CP-VRP inapplicable to the MoP-VRP, thus we need to find a new suitable algorithm to solve the MoP-VRP.
The proposed MoP-VRP considers a limited homogeneous vehicle fleet. A certain number of machines are installed in each vehicle. Production starts at the time the vehicle departs from the depot. Each customer has a strict starting time window and a soft ending time window. The vehicle should return to the depot before a certain time. The objective is to minimize the weighted sum of travel and delay cost, where the travel cost is based on the total travel length and the delay cost is based on the delayed time.
In this paper, we introduce the MoP-VRP and analyse its properties. We present mixed integer programming models for the MoP-VRP and the CP-VRP. We develop an Adaptive Large Neighbourhood Search (ALNS) algorithm for both problems. By analysing the computational results for the two problems, we illustrate the advantages of mobile production.
Most of the MoP-VRP's applications will be in the future, when 3D printing technology has advanced to a state where popular consumer goods can be 3D printed on demand. However, an application that, in principle, could be rolled out at the time of writing is that of producing spare parts needed for repairs. Most of the machinery and equipment that we rely on in everyday life contains plastic parts that can break down.
Examples of such equipment include everything from coffee machines, printers, and photocopiers to plumbing parts or even cars, buses, and planes. When a plastic part breaks down, the user may be unable to use the machine or equipment and have to wait for a repair. That repair may have to wait until spare parts are shipped from a distant warehouse, and for older equipment the spare parts may not even be available. If the spare part can be printed by a 3D printer mounted in a vehicle, the customer may get the spare part on the day she makes her request or even within hours if the MoP-VRP is solved in a dynamic setting (the present paper only considers the static MoP-VRP; the dynamic version is left for future work).
3D printers have already been used on board ships to print spare parts for repairs while at sea (Krassenstein, 2014). One may be concerned about the quality of printing as the production environment is likely to be unstable both at sea and on the road. To solve this potential risk, Phillips et al. (2020) used passive stabilization to help on board 3D printers provide equally good quality products as the traditional landbased laboratory. Based on this, we find it likely that 3D printing on board a moving van is feasible using appropriate passive stabilization. It should be noted, that the technology that we envision to be used in a spare parts delivery scenario is fused filament fabrication (FFF), which from the beginning is more robust to movement compared to the stereolithography (SLA) technology studied by Phillips et al. (2020).
The computational results in section 5.5 show that the cost for each delivery can be kept at a reasonable level if a decent amount of daily customers can be attracted to the service. The biggest challenge, in rolling out the envisioned MoP-VRP spare-parts service may be in building a vast library of 3D models that correspond to possible spare parts and in obtaining the rights to print spare parts for many different manufacturers.
The contribution of this work is five-fold: 1) We introduce the MoP-VRP, a new, non-trivial variant of the vehicle routing problem; 2) we study properties of the MoP-VRP and CP-VRP that allows us to accelerate heuristics for the two problems; 3) we design and implement Adaptive Large Neighborhood heuristics for both problems; 4) we propose testing instances for the two problems including a set of benchmark instances and a set of realistic instances in cooperation with 3D Printhuset; 5) we perform an extensive computational experiment to show advantages of the two forms of production and routing.
The rest of the paper is organized as follows: a literature review on related works is presented in section 2. The problem setting is described in section 3; based on this, two integer programming models, necessary assumptions, definitions, and properties for the defined problem are presented. We introduce the ALNS algorithm in section 4 and present computational results along with findings in section 5. Finally, concluding remarks are summarized in section 6.

Literature Review
This section presents a review of 1) the Central Production Vehicle Routing Problem and 2) mobile production technology.

The Central Production Vehicle Routing Problem
In order to effectively use resources, the optimization of production scheduling and transportation schemes in the classical MTO model have been widely studied. Different works may deal with unique problem settings, and the name given to the solved problem can vary from work to work. In the review by With respect to optimization strategy, many enterprises apply a decomposition method in practice to consider production scheduling and transportation schemes independently and optimize them separately The mathematical model and valid inequality for the CP-VRP have been widely studied in the literature. Chang and Lee (2004) used mathematical derivation to prove that the CP-VRP problem, with only one production machine, is an NP-hard problem. Geismar et al. (2008) proved that the CP-VRP problem is an NP problem and gave the lower bound of the solved problem by mathematical proof. Armstrong et al. (2008) proved that the use of the basic properties of the tackled problem can accelerate the branch-and-bound search.
Chen and Vairaktarakis (2005) defined CP-VRP with different scenarios and gave mathematical proof of the complexity of some of the defined scenarios.
The exact algorithm for CP-VRP has also been investigated in the literature. Amorim et al. (2013) verified the importance of introducing lot-sizing by solving the mixed integer programming model. Tamannaei and Rasti-Barzoki (2019) dealt with a problem with production due date and used the branch-and-bound algorithm to solve the problem. In this work, the "earliest due date first production" scheduling scheme is proposed and is proven to effectively accelerate the search for the lower bound. Dayarian  problem that jointly optimizes production, inventory, and routing decisions. The PRP usually deals with a multiple-period planning horizon, and in each period the number of products to be made, the number of products to be stored in the depot (i.e, the inventory), and the delivery plan must be decided. Different from the PRP, Both the CP-VRP and the MoP-VRP in this work will not consider overproduction (i.e., no inventory cost), and we only deal with a one-day planning horizon.

The Mobile Production Technology
To the best of our knowledge, only two works introduce the concept of mobile production technology in the literature; different works define a different kind of mobile production.
Malladi et al. (2020) introduced a mobile production mode in which production is conducted by a movable production unit called a "module". The number of modules determines the capacity of production, and the target is to optimize the allocation of modules, production, and inventory planning to serve uncertain demand within a certain planning horizon.
In Pasha et al. (2020), the mobile production is named the "factory in box". First, trucks pick up raw materials from the supplier and semi-produced products from the factory. Then, production takes place at the factory or customer's location. Although the vehicle carries the "factory" and moves during the delivery process, production still takes place at a fixed point, not en route. The mobility of production in this work is reflected in the flexible location of the factory. Production in these two works takes place at a stationary point, whereas production in the MoP-VRP takes place en route.

Problem Description and Mathematical Model
This section will describe the setting of the MoP-VRP and the CP-VRP. Then, the mathematical model for both problems will be given. After that, some properties of the MoP-VRP and CP-VRP, as well as corresponding proof are provided.

Problem Description
This section describes the setting of the MoP-VRP and the CP-VRP. We consider both problems within the context of urban delivery.
Both the MoP-VRP and CP-VRP are defined on a graph as follows: as the set of nodes and E as the set of edges defined between connected nodes. Node 0 is the depot and C = {1, . . . , n} is the set of customers. There exists a homogeneous set of vehicles K = {1, . . . , κ}, each with capacity Q. We consider a fixed number of vehicles in this work and multiple trips are not allowed. The vehicle will start from the depot and visit a sequence of customers, then go back to the depot no later than a certain time. The vehicle can depart at time 0 and should go back to the depot no later than D. Each customer i ∈ C has demand d i and production time p i . Each customer also has a pre-specified time window We consider a strict starting time window and a soft ending time window for both problems. The service time for each customer i is denoted as e i . The distance from i to j is denoted as c ij , and the traveling time from i to j is denoted as t ij . We denote M = {1, . . . , m} as the set of machines on each vehicle that produce the goods for the customers. For the CP-VRP, we define M = {1, . . . , m } as the set of machines in the depot. In this work, we assume that each vehicle has the same number of machines, so m = m * κ. In practice, there is a setup/clean up time associated with each job, but this can be included in the production time and is therefore not modelled explicitly.
For the MoP-VRP, a number of machines are installed in each vehicle. Production and delivery start simultaneously at the beginning of the planning horizon. The service start time at each location is the largest value among the arrival time, the production completion time, and the starting time window for that location.
With the current technology, there is a need for human intervention every time a new printing job is started.
This means that the driver will have to stop the vehicle and spend a minute or two starting the new job, which could be problematic while driving on the highway. This is an issue that is not modelled currently, but is left for future work. For the CP-VRP, the machines are installed in the depot, and each vehicle leaves the depot only when all the products assigned to the vehicle are ready.
The aim for both problems is to determine an optimal integrated production and delivery schedule such that the weighted sum of the total travel distance and the total delay is minimized. Clusters A and B in Figure 3 show the optimal solution for the CP-VRP and the MoP-VRP, where production starts at 0 in both problems. In the MoP-VRP, the vehicles leave the depot at 0 and production takes place on the way to the customer. However, in the CP-VRP solution, production takes place at the depot, vehicles have to wait until the production is finished before they can leave the depot. It is found that in the CP-VRP vehicle 1 will leave the depot as soon as the assigned products (1, 2) are complete, while vehicle 2 needs to wait until product 4 is made. The optimal solution for both problems has the same travel distance of 32 (5 + 6 + 5 + 5 + 6 + 5), the MoP-VRP just has a delay of 15 (2 + 6 + 1 + 6), and the CP-MoP has a delay of 54 (4 + 15 + 12 + 23). Therefore, as we set the weight for both terms to 1, the objective value for the MoP-VRP is 47 (32 + 15), which is much lower than that for the CP-VRP (86 = 32 + 54).
One way to reduce the delay to the customer could be to allow the machines at the depot to start production at an earlier time, say H time units. Therefore, we introduce early production to the CP-VRP to help reduce delay costs. Note that the vehicle should depart the depot no earlier than 0, even if all the assigned products are finished before 0. The cluster C in Figure 3 shows the optimal solution for the CP-VRP, with early production considered, in which production begins 12 time units before 0 (i.e., H = 12). It is found that the total delay can be largely reduced to 12 (0 + 1 + 5 + 6) and is lower than that in the optimal solution for the MoP-VRP.
However, the drawback of early production is that: the entire operation now spans a longer time period. We return to the consequences of this drawback in Section 5.

Mathematical Model
In this section, we present mixed integer programming (MIP) models for both the MoP-VRP and the CP-VRP.
In addition to the parameters introduced in Section 3.1, we define δ + (i) and δ − (i) as the set of nodes that can be reached from i ∈ V and the set of nodes that can reach i ∈ V , respectively. We define a dummy start node o and dummy end node d at the depot 0.
The binary variable x k ij equals 1 if vehicle k ∈ K travels from i to j, and 0 otherwise. Variable y i denotes the delay time for customer i ∈ C. The non-negative variable s k i is the time vehicle k ∈ K which starts delivery at node i ∈ V . In the MoP-VRP model, the non-negative variable v k il is the production starting time for node i ∈ V on machine l ∈ M in vehicle k (if i is served by k) and the binary variable w k ijl equals 1 if customer j's order is produced immediately after customer i's on machine l ∈ M in vehicle k ∈ K.
The mathematical model for the MoP-VRP is set as: The objective (1) is to minimize the weighted sum of travel and delay cost, where travel cost is based on the total travel length and the delay cost is based on the delayed time. In this work, we set the weight for both terms as 1. Constraints (2) ensure that each customer is visited exactly once by one vehicle.
Constraints (3) guarantee that each vehicle starts and ends at its starting and ending depots. Constraints (4) are the flow conservation constraints. Constraints (5) ensure that the capacity is not exceeded for each vehicle. Constraints (6) make sure that if a customer i is visited by vehicle k, then this order should be produced on one of the machines on k. Constraints (7)-(8) are the production flow conservation constraints for each machine. Constraints (9) determine the start production time for each order. If order j is produced immediately after order i by machine l in vehicle k (i.e., x k ij = 1), then the start production time of order j should be greater than or equal to the start production time of order i plus the production time of i, which is (10) and (11) ensure that the delivery can start only if the vehicle has travelled to the customer and the production of that order is finished. We suggest thatM should be a value no smaller (13)  In the CP-VRP, w k ijl and v k il changes to w ijl and v il because production does not take place on the vehicle. v il can have a negative value because early production is allowed in this problem.
In the CP-VRP, constraints (20)- (22) should be interpreted like (7)-(9) in the MoP-VRP. Constraints (19) ensure that each customer's demand should be produced by some machine. Constraints (23) guarantee that each vehicle can only start when all the assigned products are finished. Constraints (24) set the early production.

Properties of the MoP-VRP and CP-VRP
It is easy to see that MoP-VRP is an NP-Hard problem: we consider an instance for the MoP-VRP, where the production time for each customer is 0. It is sufficient to show that to solve this special instance of an MoP-VRP is equivalent to solving a VRP.
Next we turn to a simple property of the MoP-VRP that allows us to discard a part of the feasible solution space when designing heuristics for the problem. To explain this property, it is convenient to start with a small example for an MoP-VRP with two machines per vehicle. Consider the vehicle route [1 2 3 4 5 6], where each element is the ID of a customer. A possible machine schedule for this route is Machine 1: [1 2 5] and Machine 2: [4 3 6], where the elements of the vectors indicate which customers the machine is producing for. We propose that the schedule for machine 1 is in line with the route since the machine completes jobs in the same order as deliveries are being made to the customers. Machine 2's schedule, meanwhile, is not in line with the route because the jobs for customer 4 and 3 are produced in the opposite order of their appearance in the route.
It is obvious that changing the order of customers 3 and 4 in the schedule for machine 2 can never make the solution worse: when product 4 is produced before 3, neither of the two customers can be served before the production of both orders is finished. When product 3 is produced before 4, customer 3 can be served once product 3 has been finished while customer 4 cannot be served before both products are done; this implies that with production schedule [3 4 6], the vehicle will serve customer 3 and 4 (or the customers thereafter) at the same time or earlier than production schedule [4 3 6]. This is formalized into the following proposition: There exists an optimal solution for the MoP-VRP, such that the schedule of each machine is in line with the route of the truck carrying the machine.
Proof Suppose that s * is the unique optimal solution, where customer 2 is served after customer 1, but his product is produced before customer 1 on the same machine, and another solutionŝ where the sequence of production and delivery are the same for two customers. We define f (s) as the objective value for solution s. As the route is fixed, the service for customer 1 inŝ will start no later than that in s * , and the service starting time for customer 2 is the same in both solutions. Therefore, f (ŝ) ≤ f (s * ) ; thus, s * is not the unique optimal solution. This is contradictory to the assumption.
The following proposition states that the CP-VRP can be decomposed into a production planning problem and a routing problem when the CP-VRP only involves a single vehicle Proposition 2. When there is one vehicle and m machines in the CP-VRP, the optimal solution can be found by first minimizing the production completion time and then solving a TSP.
Proof As multiple trips are not allowed in this work and the vehicle cannot leave before the completion of all the assigned items, and since the production order does not influence the route of vehicle, the problem can be split in two as indicated in this proposition.
We next investigate the structure of the schedule for a single machine in the CP-VRP. The goal of this analysis is to allow us to rule out certain machine schedules from the set of schedules to consider. This will allow us to design faster heuristics in section 4. We note that the machine schedules influence the overall solution through the times when products are available. All products to be delivered by a single route must be available before the vehicle can depart; poor machine scheduling will therefore lead to unnecessary delay costs.
The key insight is that if the schedule for a particular machine contains two or more products that are to be delivered on the same route, then these products can be produced in direct succession. For example, consider a solution where the production schedule for a machine in the CP-VRP is [ 1 (2) where each element in the vector has the form i (r) , in which i is the ID of the customer who will receive the product and r is the route that this customer is assigned to. Starting from the back of the schedule, we can gather the two products that are to be delivered by route 2. This gives the schedule , which is at least as good as the starting schedule because the products for route 2 are finished at the same time as in the original schedule while the products for route 3 and 4 are finished earlier. We now continue backwards through the schedule to the next "scattered" products which are the products to be delivered on route 4. We gather these to obtain the schedule In this schedule, the products for route 2 and 4 are finished at the same time as in the previous schedule while the products to be delivered by route 3 are finished earlier. We now have a schedule where products that are to be delivered by the same route are in direct succession.
Lemma 1 below formalizes the process of gathering two products destined for the same route. In this Lemma, we let r(i) denote the route that serves customer i. The proof of the Lemma follows directly from the arguments above.
where θ is the length of the schedule and i j indicates the customer ID that the machine is producing for. If there exists α, β, γ such that 1 ≤ α < β < γ ≤ θ and such that r(i α ) = r(i γ ) and r(i α ) = r(i β ), then job i α can be moved to position γ − 1 and the jobs that previously occupied position α + 1 to γ − 1 will now occupy position α to γ − 2. This change will not delay departure for any routes and may reduce the departure time for some routes.
From Lemma 1 we derive the following proposition: There exists an optimal solution for the CP-VRP, where the production schedule for each machine finishes all jobs delivered by one route before switching to producing jobs that are to be delivered by a different route.
Proof Consider any optimal solution to the CP-VRP instance. If the optimal solution does not satisfy the conditions of proposition 3, then we process each machine schedule that does not satisfy the conditions by repeatedly applying Lemma 1. We go through the machine schedule from the back and gather jobs destined for the same route by repeatedly applying Lemma 1. After having applied Lemma 1 η times, we are sure that the η last jobs are grouped according to their delivery route and we therefore at most need to apply Lemma 1 as many times as there are jobs in the schedule. Applying Lemma 1 does not make the cost of the solution any worse so we will therefore eventually end up with an optimal solution in which jobs are grouped according to their delivery route for each schedule. This shows that an optimal solution with the desired properties exists.
The following proposition compares the objective value of the MoP-VRP and the CP-VRP when there is exactly one vehicle available.
Proposition 4. In the case when there is one vehicle and m machines available, the objective value of the optimal MoP-VRP solution is always better than or equal to the objective value of the optimal solution to the corresponding CP-VRP instance, assuming that the early production parameter H is set to 0.
Proof If the route and production schedule from an optimal CP-VRP solution is copied to an MoP-VRP solution, the production completion time for the customers is the same for the two problems. Since the vehicle in the CP-VRP will start from the depot after the completion of all the products, the service starting time for each customer in the MoP-VRP cannot be later that in the CP-VRP. Thus, with the same solution, the MoP-VRP will never obtain a higher total cost than the CP-VRP because the travel cost is the same and the delay cost is at least as good as that in the CP-VRP.
When the number of vehicles is greater than 1, the optimal solution to the CP-VRP can be better than the optimal solution to the MoP-VRP, even when the early production parameter H is set to 0. To illustrate,  The pseudo-code of our ALNS is presented in Alg.1. In line 2, an initial solution is constructed by a greedy insertion heuristic. In lines 6-8, we choose whether to use the noise or not based on the previous performance; after that we choose the destroy and repair operator. In the solution searching phase, we apply the threshold acceptance method proposed by Dueck and Scheuer (1990) which has been shown as a good criterion for the ALNS (Santini et al.,2018). We always accept the better solution but also accept a deteriorated solution if it is within a given threshold T of the best solution found so far. The parameter T is initialized in the beginning of the search and decreased linearly to 0 at the end of the search. We will tune the initial threshold T initial in our final test to ensure the best performance of the algorithm. The ALNS stops when a certain number of predefined iterations, N max , is reached.
To enable the comparison between MoP-VRP and CP-VRP, we also apply ALNS to the CP-VRP. The pseudo-code of Alg.1 also applies to the CP-VRP.

Construction of Solution
We used the parallel insertion heuristic to construct the initial solution. We start with a set of empty routes. In each iteration, we evaluate the cost of inserting each unassigned customer to each possible route position as well as each possible production position in the assigned route. We then identify the minimum insertion cost for each unassigned customer, select the one with the smallest minimum insertion cost, and insert it to the best position. The procedure is repeated until all the customers are included in the routes.
Different from the traditional VRP, in the MoP-VRP, not only does the routing position need to be decided for each customer but also the machine production schedule for the corresponding customer demand. For each inserted routing position, one needs to examine all the possible production schedules for each machine.
However, as proven in Proposition 1, there exists an optimal solution such that the production sequence on the same machine is consistent with the route sequence. Following Proposition 1, we only consider the production whose sequence is in line with that of route in our implementation. This helps to reduce computational time significantly. The pseudo-code for the insertion strategy is shown in Alg. 2. for each Machine l on Route r do 6: Find the last customer whose delivery position in route r is before ρ 7: Compute the delay increment in route r for the customers after position ρ 8: (i * , ρ * , l * , r * ) = (i, ρ, l, r) if the total cost increment is the least

Destroy and Repair Operators
In this work, we propose six destroy operators and four repair operators. The MoP-VRP and CP-VRP share the same destroy and repair operators introduced in this section.
Common to all destroy methods is that we first randomly pick a value Φ ∈ [ λ 1 n , λ 2 n ] as the number of customers to remove. For the interval [ λ 1 n , λ 2 n ], n is the total number of customers and λ 1 and λ 2 are the minimum and maximum ratios, respectively. We will tune the lower bound and upper bound (i.e., λ 1 and λ 2 ) in our final test.
The random and worst removals use the same strategy as that introduced by Ropke and Pisinger (2006).
The random removal picks Φ customers at random to remove. The worst removal removes Φ customers with the highest cost savings ξ i . Inspired by the worst removal, Worst-Delay Removal removes Φ customers with the highest delay cost savings, and Worst-Dist Removal removes Φ customers with the highest travel distance savings. Inspired by the existing Shaw removal (Shaw,1997;Shaw,1998), we also propose two similar destroy operators, Geo Removal and Demand Removal, which first remove a random seed customer and then Φ − 1 customers that are close to the seed customer in terms of distance and demand, respectively.
Details of the destroy operators are presented in the Appendix.
For the repair operators, the proposed ALNS applies the regret-k method. The idea is to choose the insertion that we will regret most if we do not perform it. Let r ih indicate the route where customer i has the hth lowest insertion cost and ∆ i,r ih be the inserting cost of customer i on route r ih . The regret value The regret-k method will choose to insert customer i that maximizes: arg max{g * i } and will put it into the position with the lowest cost increment. Ties are broken by selecting the request with the best insertion cost. In this work, we implement the regret-1/-2/-3/-4 method, regret-1 method is the greedy insertion introduced in section 5.1 due to the tie-breaking rule.
We also introduce noise to give randomness to the objective value every time we calculate the insertion cost of a request. We calculate the noise by applying the same strategy used by Ropke and Pisinger (2006).
Each removal and insertion heuristics is given a weight that can be adaptive adjusted to measure how well each operator performs on each instance. We also give the automatically adjusted weight to the noise so that the algorithm itself can decide whether to apply the noise or not in each iteration. A higher weight indicates a more effective operator, and an operator with a higher weight is more likely to be chosen in each iteration. The mechanism and parameters to adjust the weight are exactly the same as Ropke and Pisinger (2006).

The ALNS for the CP-VRP
The ALNS for the CP-VRP shares the same framework as the ALNS for the MoP-VRP, and the destroy/repair operators for the MoP-VRP also apply to the CP-VRP. We introduce some new strategies to accelerate the ALNS for the CP-VRP.

Construction of Solution
We also apply the parallel insertion heuristic to construct the initial solution for the CP-VRP. We start with a set of empty routes. In each iteration, we evaluate the cost of inserting each unassigned customer to all the possible route positions and fix the route position with the smallest cost increment. We then identify the smallest total cost increment by trying all the possible positions in the production sequences at the depot. We select the customer with the smallest minimum insertion cost and add it to the route and machine schedule in the solution. The procedure is repeated until all the customers are included in the solution.

The Potential Production Sequence for the CP-VRP
When we identify the best production position for a customer whose routing position is already determined, we should examine all the possible positions in the production sequence for each machine to ensure that the total cost increment is the global minimum. However, testing all the possible positions can make the computation quite time-consuming. From Proposition 3, we know that grouping production for customers who are assigned to the same vehicle to a shared machine can lead to an optimal solution. We can use this property to eliminate the cost computation of a large number of unnecessary production positions, thereby reducing the computational time significantly.
Then, we provide an example to show how we narrow the search for potential production positions.
Suppose that a new customer has already been selected to be added to route r, and we now need to consider when the corresponding production should take place for each machine. Assume the machine currently produces products for five customers. Let's use a vector to indicate the route index of each customer in the corresponding production sequence. For example, [4 4 2 2 1] means that the first two customers for whom the machine produces are delivered to by vehicle 4, the next two customers by vehicle 2, and the last one by vehicle 1. Then, we may meet two cases: The reason is that for both A and B, the ready time for routes 1, 2, 4 is the same. A, however, will lead to a higher ready time for route 3. This is because for both V and W , the ready time for route 4 is the same but V can lead to an earlier departure time for route 2, which is cost saving.

The Piece-wise Linear Function for the CP-VRP
An insertion or removal of one product in the production schedule may affect the departure time for several vehicles. If the departure time is changed, the delay cost for the route may change. Further, it is time-consuming to go through all the positions in the affected route to calculate the change in the delay cost.
When the route is fixed, we find that the total delay of the route can be expressed as a non-decreasing piece-wise linear function (Delay(ψ)) of the route departure time, where input ψ is the departure time and Delay(ψ) is the corresponding delay cost. The function has a typical shape, as given in Figure 6. With the help of Delay(ψ), we can store the information of how the delay cost changes with the departure time.
Then, each time we obtain a new departure time, we can easily obtain the new delay cost and thus the computational time can be largely reduced.
To introduce how the Delay(ψ) helps save computational time, we present an example of a route in Figure   5, and its corresponding piece-wise linear function in Figure 6. In Figure 5, the two boxes illustrate the start and end of the route, and the circles illustrate the customers visited. The numbers in square brackets above the customers are the time windows. The numbers under the edges are travel times, and the numbers in parentheses are service times. This route would not be feasible in an ordinary VRPTW setting because it would be impossible to visit customer 2 after having visited customer 1 with the given time windows.
However, in both the CP-VRP and MoP-VRP, such a solution is feasible, because delay is allowed. In Figure   6, we see that the slope is different in each interval; this represents the different number of customers whose delay costs change with the departure time. The numbers shown in the x-axis, 15

Comparison of different insertion strategy
The most accurate way to calculate the cost increment when we insert a customer is shown in Alg. 3, by which we need to try all the positions in each route. For each delivery position, we need to try all the potential production sequence to ensure that the final customer insertion causes the smallest increment cost.
We call this strategy the integrated strategy.
In line 5 in Alg. 3, we need to update the piece-wise linear function for route r because the route has changed and the pre-stored cost function cannot be used anymore. We update the function temporarily, for each Route r do 3: for each position ρ in route r do 4: calculate travel increment at position ρ

5:
Update the Piece-wise Linear Function for route r temporarily 6: for each machine l do 7: for each position ν in the schedule of l do 8: Compute the delay increment for the whole solution 9: (i * , r * , ρ * , l * , ν * ) = (i, r, ρ, l, ν) if the total cost increment is the least and it will turn back to the original form when we finish the computation for position ρ in route r. We also apply the worst-case analysis to estimate the computational complexity of Alg. 3. It is found that the computational complexity for Alg. 3 is O(n 2 ( n κ + m κ 2 log( n κ + 1))) = O( n 3 κ + n 2 m κ 2 log( n κ + 1), where n is the number of customers, m is the total number of machines in the depot, and κ is the number of vehicles in the instance. In the computational complexity, n 2 , n κ , and m come from lines 1-3, line 5, and line 6, respectively. The n κ comes from the analysis that to have the same number of customers in each route will cause the highest computational complexity. The κ 2 log( n κ + 1) is composed of two parts: κ and κlog( n κ + 1), where κ comes from line 7 because there will be at most κ positions to consider due to the analysis from To save the computational time, we apply a decomposition strategy for the CP-VRP. The pseudo-code for the proposed strategy is shown in Alg.4.
It can be seen from Alg. 4 that we fix the delivery sequence at first without considering the change of ready time. We then find the best production schedule that minimizes the cost increment for all vehicles.
Computational complexity is largely reduced because we separate the loop for the routes and the machines.
By applying the worst-case analysis, the computational complexity for Alg. 4 is O( n 3 κ + nm κ 2 log( n κ + 1)), where n 3 κ comes from lines 1-5 and nm κ 2 log( n κ + 1)) comes from lines 7-12. Under the assumption that κ = O(n), the computational complexity can be simplified to O(n 2 + mn 3 ). We can omit the O(n 2 ) part, so Algorithm 4 The decomposition strategy 1: for each non-inserted Customer i do 2: for each Route r do 3: Find insertion of Customer i in route r (assume that route ready time is not changed by insertion)

4:
Save insertion if total cost increment is the best 5: end for 6: Update the Piece-wise Linear Function for the best route 7: for each Machine l do 8: for each position ν in the schedule of Machine l do 9: Compute the delay increment for the whole solution We need to redefine the approach to find k best routes in regret-k insertion because it is different from that in the integrated method. As shown in Alg. 5, for each customer, we choose the k routes with the lowest travel increment first. Then, we update the cost increment for the picked k routes by searching for the best production schedule. After that, we sort the k routes in ascending order of the updated cost increment.

Computational Results
To evaluate the performance of the proposed mathematical model and ALNS, a series of experiments have been conducted on two different instance sets. The codes have been written in a workstation by using C++ with Windows 10, Intel Core i9-7940X, 3.10GHz, 32 GB RAM. CPLEX 12.8 is used to solve the MIP models.

Instances
We generate two sets of instances. One set is based on the well-known VRPTW benchmark Solomon Instance (Solomon, 1987) and the Gehring and Homberger's extended VRPTW Instance (Gehring and Homberger, 1999), where we use the same vehicle capacity, customer locations, duration, time windows, service times, and demands as in the benchmark instances. The other set is constructed to mimic a real-life scenario where mobile production is employed.

Algorithm 5
The best k routes for each customer in the decomposition strategy 1: for each Route r do 2: Find insertion of i in route r (assume that route ready time is not changed by insertion)

3:
Save insertion if travel increment is best 4: end for 5: Find the k routes with lowest travel increment 6: for each route r ∈ the k Best Routes do 7: for each Machine l do 8: for each position in the schedule of Machine l do 9: Compute the delay increment for the whole solution To set the maximum number of vehicles in each instance, we run a greedy sequential insertion algorithm that inserts all customers in a solution. The number of routes this algorithm terminates with is the maximum number of vehicles of the instance.
For the realistic instances, we consider a delivery scenario in Copenhagen, Denmark as the setting. We generate data using the spare-parts-for-repairs case, as described in the introduction.
As shown in Figure 7, we randomly generate 100 coordinates in a 20 by 30 km rectangle in the greater Copenhagen area. One of them is chosen as the depot and the rest as customers. The distances between the positions are generated by using Google Maps API, and the travel times are calculated using the distances divided by an average speed of 50 km/h. We assume that the drivers, working hours are from 11 am to 9 pm (i.e. 600 minutes) per day.
We also generate 25 and 50 customer instances, as the company may not receive many orders at the beginning stage of this new business mode due to immature technology or the small market. For the 25 and 50 customer scenarios, we randomly pick 25 and 50 customers, respectively, from the 99 customers. The service time for each customer is a random number within 1-5 minutes.
We create six scenarios in total, each differing in terms of the type of production time and time windows. The instances are available at https://zenodo.org/record/4892330#.YLczY4Xitdg.
During the implementation, we define the early production time as H = P mκ , where P = n i=1 p i is the total production time and is the early production coefficient. We set = 0.75 and will explain why this value is chosen in Section 5.4.

Parameter Tuning
In this section, we report results from parameter tuning on the initial threshold and the number of customers to remove, which from our experience, are very important to the algorithm. We randomly pick five instances for each size (25, 50, 100, 200) (i.e., 20 instances in total) to conduct the parameter tuning test, and we let the ALNS run 10 times for each instance for each parameter setting. The values shown from Table 1 to Table 6 represent the gap between the average objective valuez over 10 runs for each parameter and the best objective value z * found in the 10 runs among all the parameters, calculated asz −z * z * × 100%.
During the tuning for the initial threshold, the range to remove customers is set as (25%-50%). Table   1 and Table 2 show the influence of the initial threshold T initial . When T initial is too low or too high, the heuristic has a higher chance of being trapped in one suboptimal area of the search space. The results show that T initial = 10% and T initial = 17.5% are the best choices for MoP-VRP and CP-VRP respectively.  When we tuned the number of customers to remove, we first test the Solomon 25,50, and 100 instances, then pick the best five ranges to continue the tuning for 200 customer instance. We tune in this way to save computational time because it is very time-consuming to test 200 customer instances. For the MoP-VRP, as shown in Table 3, the best five ranges for the instances with up to 100 customers are (20%-50%), (10%-40%), (10%-45%), (25%-45%), and (25%-35%). So they are picked to continue the tuning. Table 4 shows the average results for these five ranges on 25, 50, 100, and 200 customer instances. We find that with the increase of λ 1 , the average gap widens; it is also not good to set λ 2 as a too large (50%) or too small (35%) value. Among these ranges, (10%-40%) obtains the smallest gap and therefore is selected as the range of removing customers for the MoP-VRP.  For the CP-VRP, based on the results shown in Table 5, we pick the best five ranges as follows: (10%-35%), (10%-40%), (10%-45%), (10%-50%), and (5%-50%) to continue tuning the 200 customer instances.

Comparison between the ALNS and the CPLEX on Small Instances
We test the MIP models for both problems on the same data set, and we give a time limit of 3600s (1h) for CPLEX to solve each instance. The instance is set with three machines in each vehicle and a production coefficient µ = 2. We generate the small-size instances by randomly picking 10, 15, and 20 customers from the Solomon 25 instances.  × 100%. As we present the results with one digit decimal, some cells with small values in the columns LB GAP and GAP appear as 0.0%.
As shown in Table 7, the proposed ALNS for the MoP-VRP finds the same or better quality solutions than CPLEX in most instances within a much shorter computational time. There is only 1 out of 168 instances where ALNS obtains a worse average result than that of CPLEX. The gap between the results by CPLEX and ALNS grows as the instance size increases. The C instances are easier to solve compared to the R and RC instances; it becomes increasingly difficult for the CPLEX to find a feasible solution within 1 hour as the instance size increases.
The proposed ALNS for the CP-VRP also outperforms CPLEX as can be seen from    Comparing the delay in the CP-VRP without early production (i.e., = 0) and in the MoP-VRP, we can see that mobile production helps reduce the total delay by more than 95%. As increases, the delay in the CP-VRP decreases dramatically. When reaches 1, all the products are nearly ready for delivery at time 0, and the corresponding CP-VRP is close to a VRPTW. When = 0.75, the delay in the CP-VRP is comparable to that in the MoP-VRP, we therefore set to 0.75 in the rest of the tests.    We run the ALNS for the MoP-VRP and the CP-VRP on the Solomon instances and Gehring and Homberger's 200 customer instances and find that the key factors that influence the objective value are the production time and the number of machines set in each vehicle. Figure 9 summarizes how the total cost changes with the number of machines and production time for the MoP-VRP. As can be seen from the figure, when the machine number is fixed, the higher production time coefficient will lead to a higher total cost; when the production time coefficient (µ) is fixed, the higher number of machines in each vehicle helps reduce the total cost. This trend becomes increasingly obvious as µ becomes larger and larger.
When the production time coefficient is 1, the influence of the number of machines on the total cost is not obvious. This shows that when the production time is small, a few machines are enough, the number of machines is not the bottleneck to obtaining a good distribution plan for the MoP-VRP. When the production time coefficient increases (i.e., 3, 5), the total cost decreases significantly, typically when the number of machines in each vehicle increases from 1 to 2. The total cost decreases by a comparatively small amount when the number of machines in each vehicle increases from 2 to 4.

Results on Realistic Instances
This section analyses the travel and delay costs for both the MoP-VRP and CP-VRP. The long term costs are then estimated, and advice on operation strategy is given.  to re-iterate that the early production coefficient was chosen to make the two problems perform similarly in terms of objective value and that the early production provides a huge advantage for the CP-VRP in terms of reducing delays (see the impact of the parameter in Table 9)

The delivery cost
The key advantage of mobile production is that it makes distribution more flexible. We find that in the CP-VRP, production has to start at least 2.5 hours (i.e., 145 mins) before the vehicles depart from the depot, and in some cases up to 6.5 (i.e. 400 mins) hours before the vehicles depart. Since the vehicles depart at 11 am, this means that almost no orders from the day of operations will be included in the CP-VRP. For the MoP-VRP, production starts when the vehicles leave the depot, which means it can accept orders until 11 am. Meanwhile, mobile production can keep an acceptable delay for customers.
Regarding the influence of the number of machines in each vehicle, the overall finding is that setting one machine on each vehicle obtains the lowest total cost in the "S W" scenario. Putting more than one machine in each vehicle lowers the total cost in most scenarios in the MoP-VRP. We only consider up to four machines in each vehicle because this is what we estimate will be feasible to fit into the vans that we used to estimate long term costs.

The long-term operation cost estimation
To put mobile production into practical use, not only should we estimate the optimal number of machines in each vehicle for good distribution but also we need to investigate how many machines in each vehicle will reduce long-term costs the most. Our long-term cost estimation is shown in Table 12, and we summarize the results for long-term costs for 25, 50, and 99 customer instances in Table 13. Considering that central production already exists in real-life practice, we only estimate the long-term cost of mobile production to check the feasibility of the new logistics mode.
The prices we estimate in this work are based on a Danish setting, and costs related to 3D printing equipment are suggested by 3D Printhuset. All the estimated costs can be found in Table 11.

Driver Printer
Hourly wage e25.  Regarding the wage for drivers, we assume that each worker works for 10 hour per day, five hours per week and 50 weeks per year. The average hourly wage is set to be e25.6, which is the hourly wage for an educated driver in Denmark plus a 30 percent addition that covers insurance, holiday pay, and indirect costs.
This means that the cost per driver a year is e63,908. If workers without taking driver education are used for the job, the wage costs will be lower. If we want to estimate the long-term costs associated with the central production mode, we need to consider that the wages for workers will be different because we need more workers in the central warehouse to handle early production; if we let the drivers handle this work, we will need to pay them more. The vehicles used in the central production mode will be a lighter type because there is no need to install printers in them, therefore the price for vehicles will be different as well.
We consider a 10-year planning horizon and the total cost for planning horizon equals the investment cost plus ten times the yearly cost. From the 10-year total cost, we can calculate the average cost per year and the average cost per order given the total number of orders per year. Table 12 shows how we calculate long-term cost. The first four columns present the scenario names,  Table 13 summarizes the proportion of each cost component in the total cost for the instances in different sizes and with different numbers of machines per vehicle. The third column presents the ratio of the total costs to that when there is only one vehicle for each data size. We can see that the drivers' salary takes up the largest part of the total cost (more than 60%) and that adding more machines to each vehicle can reduce the overall cost. Figure 10 shows how the average cost per order changes with the number of vehicles at each stage. Figure   10 illustrates that when the number of machines in each vehicle is fixed, the average cost per order gradually decreases even as the total number of customers increases. When the number of customers is fixed, the average cost per order drops dramatically when the number of machines increases from 1 to 2. Based on these findings, it is estimated that the cost for mobile production and distribution will continue to decrease as the business mode grows in popularity and technology matures.

Conclusion
This paper introduces a new variant of the vehicle routing problem, called the Mobile Production Vehicle Routing Problem (MoP-VRP), where the production takes place on the way to the customer. The MoP-VRP is highly complex, as it considers multiple products, multiple machines, soft time windows, and production schedules together within the distribution plan. For this complex problem, we propose a mixed integer programming (MIP) model as well as an adaptive large neighbourhood search (ALNS) heuristic. The computational results show that the proposed ALNS is highly efficient and can find the optimal solutions for most of the small instances within a short computational time. We also apply ALNS for the Central Production Vehicle Routing Problem (CP-VRP) and devise smart strategies to accelerate the solution process. One   the MoP-VRP does not require early production and at the same time can keep a lower delivery cost than the CP-VRP. Regarding the realistic instances, we find that the MoP-VRP is feasible in practice due to its low operations costs from a long term estimation. The cost for each delivery will gradually decrease as the customer base increases.
We propose several directions for future research. First, more constraints can be added to the basic model