Solution methods for the integrated permutation flowshop and vehicle routing problem

Article history: Received: October 15, 2021 Received in revised format: December 25, 2021 Accepted: January 29, 2022 Available online: January 29, 2022 The integration between production and distribution to minimize total elapsed time is an important issue for industries that produce products with a short lifespan. However, the literature focus on production environments with a single stage. This paper enhances the complexity of the production system of an integrated production and distribution system by considering flowshop environment decisions integrated with a vehicle routing problem decision. In this case, each order is produced in a permutation flowshop subsystem and then shipped to its destination by a capacitated vehicle, and the objective is to sequence these orders to minimize the makespan of the schedule. This paper uses two approaches to address this integrated problem: a mixedinteger formulation and an Iterated Greedy algorithm. The experimentation shows that the Iterated Greedy algorithm yields results with a 0.02% deviation from the optimal for problems with five jobs and is a viable option to be used in practical cases due to its short computational time. © 2022 Growing Science Ltd. All rights reserved.


Introduction
In recent years, researchers and practitioners have been focusing their efforts to plan decisions across different but integrated production subsystems. There is a considerable body of literature regarding the integration between production and distribution issues (Chen 2004;Moons et al., 2017). There are claims to support the need for this integration such as improved service level, usually related to the minimization of total system flowtime (Ulrich, 2013), and total system cost (Rohmer & Billaut 2015). Moreover, an integrated approach is of great interest and importance to industries that produce perishable goods (Amorim et al., 2013;Viergutz and Knust, 2014;Belo-Filho, Amorim & Almada-Lobo, 2015;Devapriya et al., 2017), chemical components (Canton, 2003;Kergosien et al., 2017) and other products with short lifespan as newspapers (Russell, Chiang et al., 2008). According to the literature, most of the applications presented so far represent manufacturing environments as a single stage production, and a little has been done regarding multi-stage manufacturing environments, such as flowshops and jobshops (Moons et al., 2017). In our research, the flowshop environment was adopted in some articles (Scholz-Reiter et al., 2011;Rohmer & Billaut, 2015;Ehm & Freitag, 2016;Ramezanian et al., 2017). These authors combine production and routing decisions to minimize supply chain related costs, such as production, inventory and tardiness.
This article presents a Mixed-Integer Problem (MIP) formulation and a constructive heuristic, both considering the two subsystems simultaneously. The Integrated Production and Distribution Problem (IPDP) studied considers a permutation flowshop manufacturing subsystem connected to a routing distribution subsystem by an unbounded inventory. The distribution is done by a single vehicle with limited capacity, which requires a pre-determined time to move between each destination. A number of clients from different locations are served by the system, and to each one of them a job is assigned. The performance criterion of the problem is the makespan of the whole sequence, which is the time elapsed between the start of the production in the flowshop and the arrival of the distribution vehicle to the origin after its last delivery.
Using the notation from Z.-L. Chen (2010), this problem is written as | | (1, ), | | . In this notation, the first term refers to the machine layout, the second to any restrictions of the production subsystem, the third to the number and capacity of vehicles and the delivery method, the fourth to the number of clients and the fifth to the performance criteria. Since both subsystems, the permutation flowshop problem (Pinedo, 2008) and the vehicle routing problem (Lenstra and Kan, 1981) are -hard, the integrated problem addressed in this article is also -hard.
To solve this problem, defined in Section 2, two approaches are adopted: a MIP is formulated in Section 3; and an Iterated Greedy (IG) algorithm is presented in Section 4. All methods were implemented and the results are shown in Section 5. Afterwards, final remarks are presented in Section 6.

Literature Review
Despite recent developments in the field, the literature brings important insights to define an IPDP. Two main groups of IPDPs can be found in the literature, according to different organizational decision levels: • A lot-sizing problem approach, usually focusing on decisions of "when" and "how much" to produce to improve a cost indicator. These are tactical level decisions, covering longer planning horizons with higher aggregation of products and operations (Maxwell, 1964;Elmaghraby, 1978;Jans & Degraeve, 2008). • A production subsystem approach, dealing with manufacturing environments such as a single machine and parallel machines (Moons et al., 2017). This approach focuses on operational level decisions of "in which sequence" to produce the orders, and most of its straightforward indicators are time-and cost-related. These decisions consider more operational details and constraints regarding the orders, machines and resources involved.
These two approaches generate entirely different problems, both in mathematical formulations and in solution methods. This article uses the second approach to the IPDP, therefore the literature review is focused on which type of production subsystem is considered, as well as how the distribution is made.
Focused on the production subsystem approach, Moons et al. (2017) presented a detailed review of the state-of-the-art literature on IPDPs, and concluded that most articles consider single-stage production subsystems, either with a single or parallel machines. Regarding the restrictions of these subsystems, such as setup times and release or availability dates, most authors do not consider them in the modelling of problems, despite their impact on the production process.
Still defining an IPDP, the distribution subsystem is usually represented by two patterns: direct shipments to a single location (Ng and Lu, 2012) or a routing-enable distribution (Gao, Qi, and Lei, 2015). According to Moons et al. (2017), the majority of articles use similar vehicles with limited capacity for distribution. The number of vehicles can be pre-determined or not: Devapriya et al. (2017) used in their model an open number of vehicles that can be hired at a cost depending on the necessity, while Cheng et al. (2017) limited the number of vehicles to one.
In more recent articles, a more focused approach on the distribution of perishable goods (Devapriya et al. 2017;Kergosien et al., 2017) and the use of meta and matheuristics to achieve high quality solutions is noticed (Devapriya et al., 2017;Darvish & Coelho 2018). Cheng et al. (2017) and Kergosien et al. (2017) proposed sequential approaches to the problems, solving the production subsystem, followed by the distribution routing. Devapriya et al. (2017) prioritized the distribution decisions, since their problem involved a single product, and Darvish and Coelho (2018) showed the dominance of the integrated approach over the sequential for a supply chain problem with location, production, inventory and direct shipment distribution.
As for IPDPs with a flowshop as the production subsystem, the articles found in the literature show different variations of the flowshop problem being addressed. Rohmer and Billaut (2015) presented methods to solve integrated problems with a permutation flowshop using a sequential approach, solving the production subsystem first and using its results as inputs to obtain the distribution decisions.  (Ehm & Freitag, 2016); and non-permutation flowshop (Ramezanian et al., 2017). In all these articles, the authors approached the problems with the objective to minimize the total operation costs of the supply chain, including inventory (from both work-in-progress and finished goods), transportation and tardiness. Therefore, the proposal of this article is novel in a sense that it considers time-related factors as its objective.
Regarding the distribution systems, the cited authors consider both single and multiple-vehicle fleets, always having the carrying capacity as a constraint. The distribution in the model from Rohmer and Billaut (2015) is done by a third-party service provider, and Ramezanian et al. (2017) addressed two distribution types, direct shipping and vehicle routing. Finally, Ehm and Freitag (2016) confirmed the advantages of using the integrated approach, showing their model gave solutions with up to 30% lower costs than those obtained by solving the production and routing problems separately.
As for the use of the IG algorithm for IPDPs, Tavares-Neto and Nagano (2018) applied it for the problem with multiple parallel machines with sequence dependent setup times and a single capacitated vehicle for distribution ( | | (1, ) | | ). The authors showed that the use of the algorithm greatly improved the initial solutions obtained through a constructive heuristic and also outperformed the adaptation of a genetic algorithm proposed for a similar problem (Ulrich 2013). Finally, Abreu, Tavares-Neto and Nagano (2021) presented a biased random key genetic algorithm with hybridization of IG for open shop with routing by a capacitated single vehicle. The proposed method got competitive results.

A Mixed-Integer Model
This section presents a MIP formulation for the | | (1, ), | | problem. The objective is to minimize the value of variable , which corresponds to the makespan of the integrated schedule of the system. The restrictions are divided into four main groups: flowshop, routes, capacity and integration.
The symbols used to define a problem instance are given in Table 1, and a solution is characterized by the decision variables presented in Table 2.

Table 1
Sets and parameters used to define a problem instance Symbol Description , Indexes used to identify processing orders, varying from 0 (dummy job) to Index used to identify flowshop machines, varying from 1 to Index used to identify a route, varying from 1 to The processing time of order on machine The size (e.g., volume) of job The total capacity of the vehicle in the same unity used for The time required for the vehicle to travel the distance between and , measured in the same unity used for

Flowshop Modelling
The first group of restrictions models the functioning of a permutation flowshop. Eq. (1) calculates a Big M for the processing times. Constraints (2) and (3) guarantee that exactly one job is allocated to each position of the sequence. Constraints (4) set the dummy job completion times to zero, Constraints (5) guarantee that the machines follow the same order of jobs, and Constraints (6) guarantee that machine will only start processing a job after it is released from − 1.

Routes Modelling
This group of constraints is responsible to model the routes which the vehicle follows to distribute the finished jobs among the clients. Since the number of needed routes is initially unknown, the upper bound of routes (one for each job) is used in order to generate the constraint groups. The constraints are programmed in such a way that the routes are only activated when needed, while the inactive have all values set to zero so they are not computed in the fitness function. Constraints (7) guarantee that each location is visited exactly once. Constraints (8) and (9) force the routes to start and end at the origin, Constraints (10) maintain the routes flow along the destinations and Constraints (11) do not allow the vehicle to have the same client as origin and destination. Finally, Constraints (12) allocate the empty routes to the end of the schedule.

Vehicle Capacity Modelling
The goal of this restriction group is to guarantee that the vehicle is not loaded over its capacity. Eq. (13) calculates a Big M for the vehicle load. Constraint (14) the load of an empty vehicle to zero. Constraints (15) calculate, considering that the jobs are loaded into the vehicle in the same order that they are delivered, the capacity of the vehicle that is occupied when each job is added to a route. Constraints (16) guarantee that the total capacity does not exceed .

Production and Distribution Integration Modelling
The fourth and last restriction group integrated both systems, defining the release and arrival times of each route and the delivery times of each job. Constraints ([r1]) and ([r2]) guarantee, respectively, that the vehicle departs for route after all its jobs are finished in the flowshop and after it arrives from the previous route. Constraints ([d1]) and ([d2]) calculate the arrival of each job according to its designated position on a route.

Makespan Function
The problem aims at minimizing , which is done by using Constraints (22). Eq. (21) sets a Big M for the makespan by summing up all processing times and distances. The value of is than set to be the latest instant when the vehicle returns from a route. Constraints (22) are activated whenever = 1, meaning that is the last delivered job of route ; than the instant when the vehicle arrives back at the origin is calculated by summing the delivery date of ( ) plus the time needed for the vehicle to return ( ). The total makespan is than obtained by taking the largest of these values, which is from the last activated route.
In this article, the proposed method of using a software to solve the MIP formulation is addressed as .

An Iterated Greedy Algorithm
In this section we present an IG algorithm for the | | (1, ), | | problem, called . uses a constructive heuristic in order to obtain an initial solution, which is shown in Subsection 4.1.

Initial Solution and Insertion Mechanism
The algorithm used to generate an initial solution for the IG is an adaptation of the / heuristic from Tavares-Neto and Nagano (2018) for this problem. The / algorithm itself is an extension of the heuristic (Nawaz, Enscore, and Ham 1983) to integrate parallel-machine scheduling with the distribution subsystem. It sets the jobs in an initial sequence and applies an insertion mechanism to sequence the jobs for both production and distribution. This heuristic, named / for future reference, uses a similar mechanism, with the difference in being applied to a flowshop production subsystem instead of a parallel machines. / consists of two main phases, and its pseudo-code is presented in Figure 1: • Phase 1 -Initial Sequencing (line 1): The jobs are ordered to form an initial solution for the next phase. In this article, five different orderings were tested for the proposed algorithm. •

Phase 2 -Insertion Mechanism (lines 5-19):
The heuristic solution is constructed by inserting each job following the initial sequence. The solution is initiated containing the first job, then, one by one, the following jobs are tested in each available position, and the one which gives the smaller fitness value is kept to the next insertion.
The extension proposed for integrating the flowshop scheduling with distribution consists in applying an insertion mechanism on the routes: for each position tested in the insertion phase, the job is inserted in all possible positions within existing routes (lines 11-15) and in a new route (line 16). Thus, this algorithm produces two sequences: a sequence related to the production schedule and a sequence related to the distribution schedule.
With a fixed production schedule, the distribution problem can be approached as a single-machine production problem with release dates, in which each job is a route and the release date for each route is the time when all jobs from it are finished in the production stage. Since the objective is to minimize the makespan, the routes are then sequenced in non-decreasing order of release dates, and the makespan is calculated. Note that for each different job insertion the release dates must be updated.
J ← set of jobs in an initial order; * , S S π π ← current and best flowshop production sequences; * , V V π π ← current and best delivery sequences; Calculate completion times of jobs in ' ; S π foreach existing rout r do Insert j in the position that adds minimum distance covered by r; Re-order the routes by release date; Update values of Fitness*, * S π and * V π with the best found so far; end Insert j in a new route and place it in the position given by its release date; Update values of Fitness*, * S π and * V π if necessary; end end Return * * , ; S V π π Fig. 1. The pseudo-code for the / heuristic

The Algorithm
is a direct extension of / , allowing better solutions to be found by using local search techniques, as shown by Ruiz and Stützle (2007). As previously mentioned, Tavares-Neto and Nagano (2018) successfully applied this approach to | | (1, ), | | , showing the efficiency of this method over an evolutionary algorithm.
is divided in three steps, and its pseudo-code is shown in Fig. 2: • Initial Solution (line 1): / is applied once to construct a initial solution; • Destruction Phase (line 6): a fixed number of jobs is removed from the current solution; • Reconstruction Phase (lines 7-11): the removed jobs are reinserted in the solution using the same insertion mechanism from / ; The Destruction and Insertion Phases are repeated in sequence times, allowing the algorithm to perform the improvements on different subsections of the current solution. In each iteration, the number of removed jobs is = × , in which is a percentage of the total of jobs in the problem. Section 5 shows how both and values are obtained. Whenever the solution found is better than the current one, it is updated with the former and used as the starting point in the following cycle. Both number of cycles and number of jobs to be removed are fixed parameters of the algorithm. * , P P π π ← current and best IPDP solution, obtained by I/OTTN; ijcyc ← number of cycles of the IG algorithm; ijjobs ← number of jobs removed in each iteration; Fitness * ← Fitness ( * P π ); foreach 0<j<igcyc do Remove igjobs randomly picked jobs from P π ; Reinsert the jobs back in P π using the insertion mechanism from Algorithm 1; if Fitness ( P π ) < Fiteness* then Fiteness* = Fitness ( * P π ); * P π = P π end end return * P π ; Fig. 2. The pseudo-code for the algorithm.

Results and Analysis
For performance testing, a set of randomly generated instances was created. Since this problem was not previously studied in the literature, no benchmark databases are available. The parameters for instance generation were taken from Tavares-Neto and Nagano (2018), with the exception of the number of machines in the flowshop subsystem, for which the range from the benchmark of Taillard (1993) was used. Five different orders were tested for the initial sequence. The objective of considering all possibilities is to analyse the impact of the initial sequence in the final solution according to the parameters of the problem instances.
• SPT Sequences the jobs in non-decreasing total processing times order; • LPT Sequences the jobs in non-increasing total processing times order; • NNDIST Sequences the jobs by applying the Nearest Neighbour Search from the origin; • SIZE Sequences the jobs in non-decreasing sizes order; • SIZEDEC Sequences the jobs in non-increasing sizes order; was modelled and solved using Python 3.5.0 and IBM ILOG CPLEX 12.7.0, while was coded in C++ using Microsoft Visual Studio 2017. Both experimentations were performed in an Intel Core i5-2500K CPU running at 3.30GHz, with 4GB RAM and Windows 7 Professional. Table 3 shows the results of according to the output. Due to the complexity of the formulation, it was only able to solve optimally the set of problems with 5 jobs within the time limit (3.600 seconds). For the problems with = 10 and almost half of the = 20 set, a feasible solution was found. No feasible solution was reached for the remainder of problems, and none of those with = 80 could be compiled. The compilation was not possible 1800

Results for instances with n = 5
The first analysis compares the results from to those obtained by , considering only the problems with optimal solutions found ( = 5). From the 1800 instances, found the optimal solution in 91% of them, with little variation regarding the initial order. As for the overall quality of the solutions, the best performing algorithms used the and sequences, showing a 0.12% Average Relative Deviation (ARD) to the optimal values. The values for all initial sequences are presented in Table 4.  Table 5 shows the average computational time in for each method to solve the problem instances, grouped by the number of machines. While the relative difference between the and is big, the former is still a viable method to be used for small size instances in practical cases. For the more complex instances ( = 20) takes an average of 4.18 to find the optimal solution, and on the worst case of all the 1800 problems, 37.14 .

Results for all instances
In this section, the results of are not considered since they were obtained for problems with = 5 only. For each instance, the minimum value obtained among the five different possibilities is taken as the best solution and the ARD is calculated relative to it. Tables 6 and 7 sort the ARD by the instances size ( and ) and distance parameters ( and ) respectively.
From Table 6, it is possible to observe that for small and medium sized instances ( ≤ 20), the initial sequence has little to no effect in the performance of . However, for large instances ( ≥ 40), consistently provides the best initial solution for the algorithm. A similar trend is observed in Table 7, in which provides the best solution for problems with ≥ 20 (the origin and destinations are placed within wider areas).
A Tukey's HSD test was used in order to verify the dominance of to the remaining orders. The test was done in the PASW Statistics 17.0 software, using the 9000 ARD values obtained for each initial sequence. The results are presented in Table 8, and show that is statistically superior to the others with 95% confidence. According to Columns 1 and 2, the ARD obtained with is significantly lower than those obtained with the other sequences, while these other four have no significant difference among them.   (2018). Aside from the production subsystem, the remaining parts of the problem are similar. In both methods, the parametrization results set the same number of jobs removed in each cycle ( ) in 40% of the total.
Regarding the initial sequence, Tavares Neto (2016) showed that the most efficient order was obtained by applying the Nearest Neighbour Search on the jobs sequence dependent setup times. This sequencing was not possible in this article since the || (1, ), | | problem does not contain setup times, therefore the is the dominant initial order. These results indicate that sequences that consider the relationships between jobs more effective than sequencing based on (non-)increasing values.
Having shown the quality of the solutions obtained with , a concern related to its practical use arises from the time required for its execution, given its high complexity. However, as can be seen on Table 9, even for the largest of instances ( = 80, = 20) the computational time of averages 12.63 seconds, which shows that this algorithm is a viable option to be used in real-life cases, providing high quality, near optimal solutions in feasible time.

Final Remarks
This article studied the IPDP consisting of permutation flowshop and delivery to multiple destinations done by a single vehicle with limited capacity. A MIP formulation and an Iterated Greedy algorithm are proposed in order to solve a set of 9000 randomly generated instances.
Due to computational constraints, the optimal solutions obtained through the MIP formulation are limited to problems with 5 jobs only. For these problems, yields results within a 0.12% range from the optimal solution, while finding the best value in 9 out of 10 problems. Regarding bigger problems, to which the optimal solutions are not available, the analysis shows the viability of using due to its short computational time (13 seconds in the most complex case) and suggests the use of the Nearest Neighbour search on distance values in order to generate an initial solution for the algorithm.
Overall, the presented methods are highly adaptable for different problems, indicating that they can be applied to more complex scenarios and used in real-world industrial cases.
As guidelines for future research, the following suggestions are made: study possibilities of improvements in the destruction and reconstruction phases of in order to improve the final results and reach higher percentages of optimal solutions found for small-size instances; apply the proposed Iterated Greedy algorithm and the MIP formulation for problems with additional features such as machine set-up times and vehicle loading/delivery times; analyze the effect of using of multiple vehicles instead of a single one in the makespan of the schedule; analyze this problem and the proposed methods while taking into consideration cost-related performance measures.