A Differential Evolution Heuristic for Integrated Production-Distribution Scheduling in Supply Chain Management

A supply chain may be considered as an integrated process in which a group of several organizations, work together. The two core optimization problems in a supply chain are production and distribution planning. In this research, we develop an integrated production-distribution (P-D) model. The problem is formulated as a mixed integer programming (MIP) model, which could then be solved using GAMS optimization software. A differential evolution (DE) algorithm is applied to solve large-sized MIP models. To the best of our knowledge, it is the first paper which applied DE algorithm to solve the integrated (P-D) planning models in supply chain management (SCM). The solutions obtained by GAMS are compared with those obtained from DE and the results show that DE is efficient in terms of computational time and the quality of solutions obtained.


Introduction
A supply chain (SC) may be considered as an integrated process in which a groups of several organizations, such as suppliers, producers, distributors and retailers, work together in order to convert raw materials into end products to be distributed to end users [1]. The two core optimization problems in a SC are production and distribution planning. In the production planning, decisions such as hiring and firing of labors, regular time and overtime production, subcontracting and machine capacity levels are made for a definite planning horizon (i.e. usually a one-year period). Whereas, in the distribution planning decisions relates to determining which facility (ies) would supply demands of which market (s).
In a traditional SC, manufacturers, wholesalers, and retailers warehouse, single or no end user P-D models. Category 5: multiple-product, multiple-plant, multiplewarehouse, multiple-end user, single-transport path P-D models.
The mathematical model of this research belongs to the second category. Therefore, the literature related to this category has been reviewed. In 1993 Pyke and Cohen [5] examined the performance characteristics of a simple integrated P-D system comprising a single work center at a factory which produce multiple products, and a single retailer. The expedite batch size algorithm is able to compute expedite and replenishment inventory control policies for the entire chain. The large number of decision variables precluded the extensive accuracy testing of the algorithm. A MIP-based P-D model was presented in for which the Lagrangian relaxation was used to accommodate the production and distribution sub-problems, and sub-gradient optimization was implemented to coordinate the information flow in a hierarchical manner. This study extends the scheme by a heuristic method which modifies distribution decisions if capacity shortages occur at the production stage. A hierarchical P-D planning approach was developed for a single multinational factory transporting multiple product families to multiple warehouses or chain stores. The approach attempts to solve the problem optimally by aggregating the time periods and product families. The obtained aggregate optimal solution for the model is disaggregated for a single period on a rolling horizon basis to reduce the problem size.
Lee and Kim developed one of the most generic P-D models in the literature and proposed a specific problemsolving procedure using a hybrid approach combining analytic and simulation methods [6,7]. In the proposed model, the first shop produces a number of parts which are used in the production of multiple products at the second shop. There is only a single plant modeled with two shops and a single stack point. The production and distribution capacity constraints in the proposed analytic model are considered as stochastic variables adjusted in accordance with general P-D characteristics obtained from a simulation model. Linear programming was used for the problem formulation minimizing the production, distribution, inventory holding, and shortage costs. An LP solver was used to implement the model and ARENA was employed as the simulation tool. This study does not investigate a multi-plant scenario and disregards a detailed production plan.
A P-D planning problem between a manufacturing location and a DC was examined by Rizk et al. [8], in which multiple parallel machine centers at the manufacturing location, economies of scale on transportation costs, as well as dynamic demand at plants and distribution centers were taken into consideration. A MIP model of the production process and three different formulations representing the general piecewise linear transportation functions were used to develop three equivalent mathematical programming models. The research investigates the impact of choosing a suitable mathematical formulation on the problem-solving time.
Nishi et al. [9] studied a distributed decision-making system for the integrated optimization of production scheduling and distribution planning. An integrated optimization model was formulated using MILP. The developed model was then decomposed into production scheduling and warehouse planning sub-problems using an augmented Lagrangian approach. A distributed optimization system was developed to solve the sub-problems by gradually generating the feasible solutions through repeatedly exchanging data obtained at each sub-system (data update) to accommodate the modifications caused by unforeseen changes. The study formulates the production processes at a single plant and considers a single warehouse comprising a number of storage areas. There are no end-users involved in this model and production/transportation alternatives are simplified.
Farahani and Elahipanah [10] solved a MILP bi-objective model for a just-in-time (JIT) distribution planning of a three-echelon SC using multi-objective Genetic Algorithms (GAs).
Two functions were considered for optimization: cost minimization as well as minimization of the sum of backorders and surpluses of products in all periods. In fact, the second objective function represents the JIT delivery and minimizes the earliness and tardiness of product deliveries. Delivery lead-times and capacity constraints were considered for a multi-period, multi-product and multi-channel network. Since the study primarily concerns the distribution of items from suppliers to retailers, this model replaces the fixed and variable production costs by the purchasing costs and completely disregards the production issues and the multiplicity of manufacturing plants.
A new solution approach was proposed by Safaei based on the integration of mathematical and simulation techniques to solve an integrated multi-product, multi-period, multi-site P-D planning problem [11]. A MILP model was developed for formulating the problem. In this study, to consider the stochastic factors, a hybrid mathematical-simulation approach was proposed consisting of independent mathematical and simulation models. The proposed model can be treated as a single-plant model with multiple machine centers and many characteristics of a real world P-D planning problem are disregarded in this model (e.g. detailed production alternatives, production cost elements, backlogging costs, and inventory management issues).
In this study, we develop an integrated productiondistribution (P-D) model. The problem is formulated as a mixed integer programming (MIP) model, which could then be solved using GAMS optimization software. A differential evolution (DE) algorithm is applied to solve large-sized MIP models. The solutions obtained by GAMS are compared with those obtained from DE and the results show that DE is efficient in terms of computational time and the quality of solutions obtained. [12].

Figure. 1. A two-dimensional cost function showing its countour lines and the process for generating
This paper is organized as follows. Section 2 explains the characteristics of the problem; in section 3, the mathematical model is developed and section 4 proposes a differential evolution (DE) algorithm to solve the model; the computational results and a discussion are included in section 5. Finally, the conclusions and proposals for future research are presented in section 6.

The Problem Definition
In this section, we explain the problem definition and basic assumptions to better understand the problem. As mentioned before, in this research we consider a multiple-product, single plant production-distribution model in (SCM). In the scheduling system considered in the production line, is flow shop. Batch processing is also applied which means that after the execution of a series of jobs (products) in the production line, the relevant jobs (products) are delivered to the vehicles center and vehicles distribute them to their destinations (customers). Customers are geographically dispersed around the depot (plant). This problem is called vehicle routing problem (VRP). Vehicles have capacity constraints and the number of vehicles are defined (VRP with capacity constraint). Each job (product) has a specific due date. We regard inventory cost when the jobs (products) deliver to customer earlier than due date and we have penalty cost at the time of tardiness.
The considered assumptions in the production part are as follows: The processing time of each job on each machine is definite. The number of jobs and production machines are deterministic.
Jobs are uninterruptable which means that the production of a product must be continued up to the end of the process. Processes are uninterruptable which means that no machine can process more than a job at the same time.
No job can process on more than a machine simultaneously.
No machine's failure exists and they are always available during the scheduling. The considered assumptions in the distribution part are as follows: Each customer must be serviced only once and only by a vehicle. Vehicles must start and terminate their rout from depot (production plant). Total demand of customers must be less than the capacity of vehicle. The last customer of each route must be serviced before the specific considered time.

The Mathematical Model
The following notations are used for the mathematical formulation of the proposed model: , , , , , Objective function (1) minimizes the transportation cost, as well as the earliness and tardiness cost of product delivery. Constraint (2) represents that one job is assigned to each situation. Constraint (3) represents that one situation is assigned to each job. Constraint (4) assures that job 1 is processed on all machines without any delay. Constraints (5-6) guarantee that there is no idle time on machine 1. Constraint (7) represents that minimization time of each job on machine l+1, won't be earlier than the completion time of job on machine l. Constraint (8) expresses that the job which is in the situation p+1 won't initialize on machine l till the job on situation p is finished. Constraint (9) shows the completion time of batch i. Constraint (10) assures that depot (production plant) is available in zero time. Constraint (11) expresses the available time of vehicle k. This is the time that all the jobs of the batch which is assigned to vehicle k, is completed. Constraint (12) represents the capacity satisfaction of vehicle. Constraint (13) depicts that when a vehicle is entered to a node, it must be sent out from the mentioned node. Constraint (14) delineates that each vehicle must be serviced only once. Constraints (15-16) express that number of vehicle for service to customer is exactly equal to v. Constraints (17-18) represent the delivery time of job i. Constraints (19-20) express the tardiness and earliness in delivery of products to customers respectively. Constraints (21-23) represent the range of model variables.
To solve the model with GAMS optimization software, it is necessary to linearize the non-linear constraints (constraints 9-17-18 are non-linear).

DE Algorithm
Differential Evolution (DE) is a stochastic, populationbased optimization algorithm introduced by Storn and Price in 1996 [12]. Global optimization is necessary in fields such as engineering, statistics and finance. But many practical problems have objective functions that are non-differentiable, non-continuous, non-linear, multi-dimensional or have many local minima, constraints or uncertainty. Such problems are difficult to solve. DE can be use d to find approximate solutions to such problems. DE uses NP D-dimensional vectors parameter similar to Eq.(27) as a population for each generation G. , , 1, 2,....., NP is the number of population for each generation which does not change during the minimization process. The initial vector population is chosen randomly and should cover the entire parameter space. We consider a uniform probability distribution for all random decisions. In case a preliminary solution is available, the initial population might be generated by adding normally distributed random deviations to the nominal solution ,0 nom x . DE generates new parameter vectors by adding the weighted difference between two population vectors to a third vector. This operation is called "mutation". The mutated vector's parameters are then mixed with the parameters of another predetermined vector, the target vector, to yield the so-called trial vector. Parameter mixing is often referred to as "crossover". If the trial vector yields a lower cost function value than the target vector, the trial vector replaces the target vector in the following generation. This operation is called "selection". In the followings, the different parts of the proposed DE for the proposed P-D model are described.
In (2) { } the amplification of the differential variation 2 3 , , Figure 1 shows a two-dimensional example that illustrates the different vectors which play a part in the generation of

Cross-over
Cross-over is introduced in order to rise the diversity of the perturbed parameter vectors. The trial vector is: ,

Selection
The greedy criterion is used in order to decide whether a vector it should become a member of the next generation, i.e. G + 1. So, the comparison between the trial vector , 1

Computational Results
This section explains the test problems which aim at showing the applicability of the proposed algorithm. The computational results are reported, evaluated and analyzed with respect to the proposed model. For the small and medium sized problems, the solutions presented by DE are compared with the results obtained from GAMS optimization software.

Designing the Test Problems
Various test problems, with different sizes are considered to evaluate the performance of the proposed algorithm. We consider two sets of 12 small and medium sized and 6 large sized problems to be solved using the DE, i.e., a total of 18 instances were run.
In each problem, the values of each group of parameters are generated randomly between their lower and upper bounds, based on Table 1.

Setting the DE Parameters
Parameters of the designed DE include CR, F, population, and iteration. Primary tests are carried out in order to determine the values of these parameters. A trade-off between the solution time and the quality determined the appropriate value. The values for CR, F, population, and iteration are set to 0.5, 2, 100 and 200 respectively.

Computational Results
The proposed P-D model has been solved using CPLEX solver of GAMS 23.6 for small and medium sized problem. For large sized problem the proposed DE is employed. The proposed DE is coded in Matlab 7.0.4. All the test problems are solved on a Pentium 4 computer with 448 MB RAM and 2.80 GHz CPU and the results are summarized in Table 2.
A quality criterion, ERROR, is defined to show the deviation of the value of the DE solutions from the values of GAMS, according to Eq.(31). .
The results are observed in Table 2. To investigate the performance of GAMS optimization software and DE heuristic algorithm, two criteria, i.e. objective value and run time for small and medium sized problems have been considered. On average, DE algorithm achieved 91.7% of the exact optimal solutions within 16.7% of the exact run time. In other words, the average deviation from the optimum for the small and medium sized problems does not exceed 8.3% of error. The trivial deviation shows the efficiency of the proposed algorithm.
We considered 1000s as a limit to stop the algorithm after a specific run time. As it is clear in Table 2, an error message is appeared in the test problems 11 and 12 run by GAMS optimization software, which means that the run time exceeded 1000s. Figure 2 depicts, DE run time versus GAMS run time. It is shown that GAMS run time increases exponentially as the problem size increases. On the average, DE algorithm achieved approximately 92% of the exact optimality in 16% of the exact run time. Figure 3 delineate DE objective value versus GAMS objective value. It is clear that DE objective value compared to that of GAMS decreases when the problem size increases. Figure 4 delineates the trend of optimality ratio and run time ratio. It is observed that decreasing in slope of run time ratio is more than the decreasing in slope of optimality ratio. In other words, the ratio of DE computational time compared to that of GAMS decreases when the problem size increases. This means that DE is efficient in terms of computational time and the quality of solutions obtained.
Objective value and run time for GAMS optimization software and DE heuristic algorithm for large sized problem are calculated. As shown in Table 3, GAMS error message for exceeding the resource limit appear in test problems 13 and 14.
Furthermore, out of memory error message of GAMS appears regarding test problems 15, 16, 17 and 18. This is because of exponential increase in GAMS run time which leads to the saturation of RAM and stopping the calculation before achieving an acceptable solution.

Conclusions and Future Research Directions
A new model for the integration of production and distribution scheduling in supply chain management has been developed in this paper which minimizes the costs of transportation, earliness and tardiness of delivery a product. We also considered batch processing in production line, which leads to cost and time minimization of the whole supply chain.
A DE algorithm is designed to solve the model for the large-sized problem instances of this mixed-integer programming problem. To the best of our knowledge, this paper is the first paper which used DE heuristic algorithm to solve the P-D planning models in supply chain management. Some small-sized and medium-sized test problems have been solved. The results obtain from DE satisfactorily compared with the results obtained from GAMS optimization software. Some other large-sized test problems have also been solved using DE.
The following approaches are proposed for future research: In this research, plant is regarded as a single-depot for distribution. Considering multiple-depots can be the subject for future research. Cost minimization has been recognized as the most respected performance measure for the evaluation of SC performance. Maximizing the service level and also maximizing the profit can be considered as P-D objective functions. Products can have expired dates. This influences the order time and the maximum duration allowed for holding the products. In this research, we considered identical vehicles for distribution of products. Consideration of non-identical vehicles with different speeds can be the future point of view.
In this research, we considered a flow shop scheduling system in production line. Consideration of job shop scheduling system can be more challenging for the future.
In this research, we only focused on the integration of production and distribution. It is also possible to add the supply-side and integrate the whole supply chain. In this research, the demand certainty is considered.
Considering the demand uncertainty can improve the applicability of the considered model. There is still a need to further extend the effectiveness of the existing solution approaches and to test the new arrivals such as Ant Colony Optimization (ACO) and Bee Colony Optimization (BCO) techniques.