Solving a fuzzy multi objective model of a production–distribution system using meta-heuristic based approaches

This paper studies a multi-objective production–distribution system. The objectives are to minimize total costs and maximize the reliability of transportations system. Each transportation system is assumed to be of unique reliability. In the real world, some parameters may be of vagueness; therefore, some tools such as fuzzy logic is applied to tackle with. The problem is formulated using a mixed integer programming model. Commercial software can optimally solve small sized instances. We propose two novel HEURISTICS called ranking genetic algorithm (RGA) and concessive variable neighborhood search (CVNS) in order to solve the large sized instances. RGA utilizes various crossover operators and compares their performances so that better crossover operators are used during the solution process. CVNS applies several neighborhood search structures with a novel learning procedure. The heuristics can recognize which neighborhood structure performs well and applies those more than the others. The results indicated that RGA is of higher performance.


Introduction
Nowadays, most firms are using optimization models, algorithms and decision support systems to improve their operations, in terms of flexibility, costs, quality of products, and distribution systems, to remain in competitive markets (Yilmaz and Catay 2006).Supply Chain Management (SCM) is a suitable and necessary tool to prosper in this competitive market.The most important purpose of SCM is integration of decisions in the different parts of a system, e.g.outsourcing, production planning, inventory and distribution management (Simchi-Levi et al. 2004).
The integration of production and distribution decisions is a key problem in supply chain networks (Gen and Syarif 2005).To integrate such decisions, we consider a threeechelon system including production, distribution and customer elements.This problem includes transporting products from production facilities to distribution centers (DC) and transporting them from DCs to customers.
There are different reasons such as high transportation fixed costs, the length of the set-up times for manufacturing the products and the constraint of the production capacity, which cause customers' demands to be sometimes supplied with delay.In these situations, incurring stock out maybe authorized; however, this assumption is not considered by most of the previous researches.Another realistic assumption, commonly ignored, is the reliability of each type of transportation system (TS) used to shift products from plants to DCs and from the DCs to customers.Few researches mention that an unreliable TS may result in a large quantity of losses while transporting products in the supply chain (Pokharel 2008).
The major objective of this research is presentation a fuzzy multi objective production-distribution (PD) model with multiple products, multiple periods, while having mul-tiple plants and multiple distribution centers.In this model, customers' demands are not fully satisfied.We classify the customers based on their unit shortage costs in such a way that the high priority is given to customers with the highest unit shortage costs.We also consider several types of TSs where each pair of TS and route (as a system) has its own reliability.Customers' demands, available time for each plant and reliability of each TS are considered as fuzzy numbers.We apply classic GA, a heuristic method based on GA called ranking genetic algorithm (RGA), classic VNS and a heuristic method based on VNS called concessive variable neighborhood search (CVNS) in order to solve the problem.
The remainder of the paper is organized as follows: in section two "Literature review" on PD system is given.Section three "Problem definition" explains the problem statement and formulation of the model.In section "The proposed solution algorithms", we propose the heuristic algorithms to solve the model.Section five evaluates the "Computational results".Finally, the last section gives "Conclusions and ideas for further researches".

Literature review
In this section, we give some selected researches on PD systems together with the solution methodologies.Sarmiento and Nagi (1999) presented a review on integration analysis of PD systems.They discussed some researches, which explicitly considered TS and other logistical aspects.Lejeune (2006) proposed a periodic three-stage supply chain network for a PD inventory system.He formulated this system by a mixed integer linear programming (MILP) model and developed a meta-heuristic algorithm based on VNS.Aliev et al. (2007) presented a multi-product and multiperiod fuzzy integrated PD model in which some parameters such as maximum capacity of DCs, maximum production capacity and demands for all products were assumed to be fuzzy numbers.Selim et al. (2008) considered an integrated PD model focusing on integration of different objectives of producers.The major objective is to arrive at a suitable profit for all producers.They designed a multi-objective linear programming model for that problem.Liang (2008) proposed a fuzzy multi-objective model with the piecewise linear membership function for an integrated multi-product and multi-period PD system.He assumed that decision makers were familiar with estimating optimistic, pessimistic and most likely parameters of the model, such as customers' demands and available machine capacities.Bachlaus and Pandey (2008) considered a multi echelon supply chain network focusing on agility as an important criterion.The network includes suppliers, plants, distribution centers, cross docks and customers.They designed a multi objective model which minimizes the total cost and maximizes the plant flexibility.A novel algorithm called Hybrid Taguchi Particle Swarm Optimization (HTPSO) is proposed to solve the problem.The results confirm the effectiveness of proposed methodology in solving the problem.Kazemi et al. (2009) have proposed a PD system including supplier, producer, assembler, distributor and customer.They applied two scenarios to solve this problem.The first scenario is a centralized GA and the second one is an agent-based GA.Jawahar and Balaji (2009) considered a two-stage distribution problem.Two types of costs were involved in the problem.The first one was continuous cost linearly dependent on the quantity transported from a source to a destination; the second one was a fixed charge, incurred whenever there was a transportation of a non-zero quantity from a source to a destination.The objective is to minimize the total distribution costs.Jolai et al. (2010) proposed a multi-producer, multidistributor and multi-product PD problem.They designed a fuzzy multi-objective model for this system and applied GA, particle swarm optimization (PSO) and a heuristic algorithm based on GA to solve this problem.Delavar et al. (2010) considered a production-distribution problem with air transportation to optimize customer service at minimum total cost.In order to solve the addressed problem, two GA approaches were developed.To illustrate the effectiveness of the proposed algorithm, different sized problems were designed and the computational results were compared together.Kuo and Han (2011) proposed a two-echelon linear programming model for a distribution problem.They applied an efficient method based on GA and PSO in order to solve the problem.
Che (2012) presented a decision methodology for the multi-echelon unbalanced production distribution planning.Production loss, transportation loss, quantity discount, production capacity, and starting-operation quantity are considered simultaneously in the proposed multi-product model.A modified PSO was proposed to solve the model.Finally, two simulated supply chain structures were proposed to study the effectiveness of the solution method.Amorim et al. (2012) presented a multi objective programming model for an integrated PD system.Products have a short and fixed life cycle and they are corruptible in this system.Bashiri et al. (2012) proposed a mathematical model for a strategic programming in a multi echelon system with multi producers and distributors.They solved several numerical examples to show the efficiency of the proposed model.Liu and Chen (2012) presented an integrated model for inventory routing and scheduling problem.Since the search for obtaining the optimum solution in this problem is very time-consuming, they used a VNS algorithm in order to solve the problem.Rajeshwar et al. (2012) considered a three echelon PD system in a supply chain network.They assumed uncertainty in demand, which in turn, results in increasing costs.This model was solved by PSO algorithm.Raa et al. (2013) considered an aggregate productiondistribution problem for a manufacturer of plastic products with injection moulding.Their proposed model is of MILP formulations.A meta-heuristic solution approach is used for solving the model.The main purpose of the approach is to quantify the opportunities which mould sharing offers to the plastics manufacturer.The results of computational experiments reveal that mould sharing can reduce the productiondistribution total cost around 10 %.Hiremath et al. (2013) considered a distribution network in a manufacturing supply chain.The network includes a set of customers with known location and a set of manufacturing plants and distribution centers with unknown location.A single type of product is manufactured by the plants.They proposed a MILP model minimizing the total cost, maximizing the fill rates and maximizing the resource utilization of the facilities.The model is solved by a Multi Objective Genetic Algorithm (MOGA).They showed the effectiveness of the applied algorithm analyzing the computational experiments.Fujita et al. (2013) studied a supply chain network with global product family.They proposed a MILP model to obtain plants locations for module production, assembly and final distribution.The objective of the model is to minimize the total cost.Bandyopadhyay and Bhattacharya (2013) considered a two-echelon supply chain network.They formulated the problem as a multi-objective programming model minimizing the total cost and the bullwhip effect.They applied the non-dominated sorting genetic algorithm II (NS-GA II) to solve the problem with new crossover and mutation operators.They compared the results with those from the original NS-GA II.Cui (2014) considered the challenging joint optimization problem of production planning and supplier selection.This integrated problem has been formulated as a MILP model.The objective of the model is to maximize the total profit of the manufacturer.An improved GA was proposed to solve the given problem.Lee and Kim (2014) presented an integrated PD model to obtain the optimal policy in a single vendor and single buyer network.The products are assumed deteriorating and defective.The objective of the model is to maximize the supply chain profit.Numerical examples and sensitivity analysis are provided to study the validity of the given model.
The fuzzy set theory proposed by Zadeh (1965) is an effective tool in order to tackle with the uncertain data.Wang and Fang (2001) and Jamalnia and Soukhakian (2009) proposed models and solution approaches for fuzzy integrated production models.Kumar et al. (2014), presented a twoechelon integrated procurement production model for a single manufacturer and a single buyer supply chain.Customers' demands are represented with fuzzy numbers.They determined an optimal multi-ordering policy for the buyer in one production cycle of the manufacturer.They proposed an algorithm to extract the optimal policies for both parties.Finally,

Distributors Customers
Fig. 1 Three echelon supply chain network they solved a numerical example to illustrate the effectiveness of the algorithm.Bai and Liu (2014) presented a new robust optimization method for supply chain network design problem.They selected a fuzzy value-at-risk optimization model with uncertain demands and transportation costs.They decomposed the model to two equivalent mixed integer parametric programming sub-models.They presented an applicable example and showed that the proposed method could provide an effective and flexible way to design a supply chain network.

Problem definition
This paper considers a three-echelon supply chain system including multiple plants, multiple DCs and multiple customers.A schema of the addressed system is presented in Fig. 1.Each plant produces a number of different products and moves them to the DCs with constrained capacities.Unsatisfied demand of each customer is backordered; however, the whole customers' demands of different products must be satisfied by the last period.We consider several types of TSs with a fixed reliability rate.We also consider a reliability rate for each route with respect to the number of the accident records, geographical conditions and so on.Reliability rate of each TS in each route is obtained from multiplication of the addressed reliability rates.
The additional assumptions are that we can only produce one batch of each product in each plant during a single period of time; each transportation system can move several times at each time period in a given route; each transportation system has a limited capacity.
After wards, we introduce the parameters, indices and the decision variables used in the model.Thereafter, the mathematical model in the form of a MILP is presented.

Parameters and indices
To formulate the problem, some distinctive parameters are required.We need to define some indices for the major elements of the model.We give six indices.

I
The number of plants

J
The number of DCs

K
The number of customers

P
The number of products

T
The number of time periods

M
The number of TSs The mixed integer linear programming model The studied problem is formulated in a mixed integer linear programming model with two objective functions.
The first objective function in Eq. ( 1) is to minimize the total operations costs, including production costs in plants, transportation costs from plants to DCs and from DCs to cus-tomers, inventory holding costs at DCs and inventory shortage costs at customers' sites.The second objective function in Eq. ( 2) is to maximize the reliability rate of transporting the products from plants to DCs and also from DCs to customers.Constraint (3) ensures that the total processing time of products is not more than the maximum available time.Constraint (4) indicates that all products manufactured at each plant must be transported to DCs at the same time period.Constraint (5) represents the inventory balance for the DCs.Constraint ( 6) is concerned with the capacity of each DC at the end of each time period.Constraint (7) controls that each plant manufactures each product at each time period if the corresponding binary variable takes value of one.Constraints ( 8) and ( 9) are concerned with the TS capacity.Constraint (10) represents the inventory shortage balance in customers' sites while Constraint ( 11) is the same as (10) but for the last period.Constraints (12-14) define the non negativity, integer and binary status of the decision variables, respectively.

Fuzzy parameters of the model
In real world situations, some parameters maybe associated with vagueness (Sakalli 2010).In the proposed model, we apply fuzzy set theory to cover this problem in order to attribute values to the imprecise parameters; the considered parameters are maximum available time for each plant (Maxt it ), customers' demands (dem pkt ) and Coefficients of reliability function (R AI m i j , R AJ m jk ).

Maximum available time (Maxt it )
We add the overtime ( t it ) to Maxt it in order to increase the flexibility of the basic time, so that if we use the overtime, we have to pay higher cost to produce the products.Constraint (3) is replaced with the objective function in (15) and the constraint given by ( 16).
where, 0 ≤ α ≤ 1. Let's consider Maxt it as the fuzzy number.We give a numerical example to indicate how to apply the given replacement.Suppose that Maxt it = (500, 560), Maxt it = 500 and t it =60; the result is as follows: Customers' demands (dem pkt ) We use triangular fuzzy number in order to tackle the vagueness of customers' demands.We apply Centroid method [as given in Eq. ( 17)] to represent the value of demand by an exact number as in Chen and Pham (2001).
where, μ(x) represents the membership function of fuzzy number x.
We consider dem 111 as the demand fuzzy numbers with the membership function μ (dem 111 ).Suppose that the fuzzy number is as (40,43,50).μ (dem 111 ) and the Centroid value are calculated as follows: R AI m i j and R AJ m jk are the reliabilities of each transportation system in each route.It is clear that considering an exact value for these parameters is not easy.We assume they can be represented by triangular fuzzy numbers.The second objective function of the model in Eq. ( 2) is converted to what is given in Eq. ( 18).
where l R AI m i j , m R AI m i j , h R AI m i j are the lower, median and upper bounds of parameter R AI m i j , while l R AJ m jk , m R AJ m jk , h R A J m jk are lower, median and upper bound of parameter R AJ m jk , respectively.To convert this objective function to crisp form, we use the given technique by Zimmerman (1983).According to this method, the distance between the lower and median bounds is minimized while the distance between the upper and median bounds is maximized.Eq. ( 18) is replaced with Eqs. ( 19)-( 21)

The proposed solution algorithms
The problem under study is NP-hard, accordingly, the exact methods cannot solve the problem efficiently.Therefore, we apply meta-heuristic algorithms to solve the problem.Meta-heuristics are approximate solution algorithms obtaining optimal or near optimal solutions within a reasonable amount of computational time (Haupt and Haupt 2004).

LP-metric method
There are different techniques such as LP-metric, maxmin, weighting and lexicography in order to solve a multiobjective problem and to convert it into a single-objective one.The objectives of the problem under consideration are contradictory, since, increasing the transportation reliability causes the total cost to increase.LP-metric is one of the well-known techniques to simultaneously consider conflicting objectives (Mazdeh et al. 2010;Mirzapour et al. 2011).This technique concentrates on minimization of deviations from the best values.To apply this technique, the ideal value of each objective function should be obtained, thereafter Eq. ( 22) is Minimized In which f * j and γ j represent the ideal value and the weight of objective function j. p is a control parameter which must take an integer value equal to or greater than 1 (Opricovic and Tzeng 2004).Assuming p = ∞, the problem turns into the minimization of the maximal deviation as in Eq. ( 23).

Penalty function
Since we will use meta-heuristic algorithms in order to solve the problem, we need to control the infeasible solutions.A penalty function is designed so that it takes zero for feasible solutions and a large positive value for infeasible solutions.
The penalty function of a solution x, called P(x), for constraint g(x) ≤ b is calculated using Eq. ( 26) (Yeniay and Ankare 2005) In which R is a large positive number.The penalty value is added to the objective function value.

Ranking genetic algorithm
GA is an optimization and search technique based on the principles of genetics and natural selection.The method was developed by Holland (1975).Some instances of GA for solving supply chain problems are given in the rest.Kannan et al. (2010) compared GA with the solutions obtained by GAMS optimization software.The solution reveals that it performs very well in terms of both quality of solutions and computational time.Sourirajan et al. (2009) indicated that GA yields significantly better solutions than a greedy heuristic.Aliev et al. (2007) employed a fuzzy GA method to solve a production distribution problem and the experimental results showed its high efficiency.
The proposed ranking genetic algorithm (RGA) in this paper can be described as follows: There is a set of solutions (or chromosome) called population.RGA iteratively evolves the population and generates a new set of individ-uals.At each iteration which is called generation, some of the best individuals are directly copied into the next population in order to reproduce the population.Then, individuals in the current population are modified by a set of genetic operators.One partially random mechanism, called selection procedure, selects two individuals from the current population.Another operator, called crossover, combines the two selected individuals and generates two newer individuals.A proportion of individuals in the next population is produced by the crossover.Another operator, called mutation, produces the remaining individuals of the next population.Finally, the current population is replaced with the individuals of the new population, and the algorithm restarts.
We use the elite strategy to protect best solutions of the former populations.Using the elite rate, we specify number of chromosomes directly copied to the next population as in Eq. ( 27) In which N R and P rep represent the number of chromosomes copied by elite operator and elite rate, respectively.To apply the crossover operator and generate new chromosomes, we need a selection mechanism which randomly selects from chromosomes of the current population.It should work so as to give more chance of being selected to chromosomes of better fitness.We use the roulette wheel mechanism for the selection procedure.
In the proposed RGA, we use several types of crossovers including one-point crossover (Fig. 2a) and multi-point crossover (Fig. 2b).Each crossover is given specified rank with respect to the quality of its solutions in each population; as the best crossover takes the highest weight and the worst crossover takes the lowest weight.The weight of each crossover can be different in various iterations, but all types of crossover operators are applied to produce new chromosomes in all iterations.The weight of each crossover shows the proportion of individuals produced by that crossover.Therefore, the best operator in previous population makes most of the new chromosomes in current population [Eq.( 28)-( 29)].
N C i is the number of chromosomes produced by ith crossover, N pop is the population size, P cr is the crossover rate and w i is the weight of ith crossover.
After applying the elite and crossover operators, the remaining chromosomes are produced by mutation operator as follows: one chromosome is selected using the selection mechanism, previously used to choose the parents for crossing.We make a slight change in the selected chromosome and generate a new chromosome using the mutation operator.The mutation operator might produce a worse chromosome.The proposed algorithm prevents to make worse solutions when applying the mutation operator.We first produce a new chromosome and then calculate its fitness value.If this value is better than the average fitness of the previous population, we copy the new produced chromosome to the next population; otherwise, the chromosome is omitted.We specify the number of chromosomes by the mutation operator by Eq. ( 30)-( 31) where, N M and P mu represent the number of obtained chromosomes using mutation operator and mutation rate, respectively.
Figure 3 illustrates the flowchart of the proposed RGA.To emphasize the difference between RGA and GA methods, the specific features of RGA are highlighted with gray color.Notice that N(repro), N(cross i ) and N(mut) in Fig. 3 are counter of obtained solutions using elite, ith crossover and mutation operators, respectively.
The concessive variable neighborhood search VNS is a meta-heuristic algorithm based on the changes of neighborhood structure.This algorithm has been proposed by Mladenovic and Hansen (1997).Generally, the local search based meta-heuristic algorithms, such as simulated annealing and tabu search, repeatedly use a constant neighborhood search structure which is the main mechanism of reduction the diversity of solutions in these algorithms; but, VNS applies several neighborhood search structures, and intelli-   The problem under study is NP-hard and exact methods can only solve small sized instances.The small sized instances can be solved by LINGO 8.0 within 100-500 s, while for the large sized instances, LINGO acquires no optimum or near optimum solution after 1000 seconds.Metaheuristic algorithms (GA, RGA, VNS and CVNS) are applied to solve all sized instances.For this purpose, we tune the parameters of GA and RGA for all small and large sized instances using response surface methodology (RSM).RSM is a mixture of mathematical and statistical methods applied to determine the best value of each algorithm's parameter (Montgomery 2004).The optimal values of GA and RGA parameters for all instances are given as in Table 3.In this table, N pop , P cr and P mu represent population size, Crossover rate and mutation rate of GA and RGA.The last column is this table gives the time limits given to algorithms to proceed.All algorithms converge within the time limits.Figure 6 gives the corresponding convergence charts.
We apply relative percentage deviation (RPD), to compare the effectiveness of GA, RGA, VNS and CVNS.To calculate RPD values for an objective function of minimization type, Eq. ( 32) is used (Naderi et al. 2009).

R P D
where Alg and Min represent the objective function value of each problem instance and the minimal objective function value of the corresponding problem.
For the maximization type, Eq. ( 33) is used.
where, Max represents the maximal objective function value of the corresponding problem.Table 4 gives the average of RPD of the problems for each algorithm with respect to the first and second objective functions of the model (Z1 and Z2).LINGO can only solve small sized problems.As is clear from Table 4, RGA and CVNS outperform the classical GA and VNS for almost all problems.
Figure 7 gives the plots of the average and interval RPD values of GA, RGA, VNS and CVNS for small and large sized problems.LINGO and RGA outperform the other algorithms for the small sized problems (as Fig. 7a) while RGA outperforms the other algorithms for large sized problems (Fig. 7b).Furthermore, RGA and CVNS outperform the classical corresponding forms (i.e.GA and VNS) for small and large sized problems.
The results are analyzed using means of analysis of variance (ANOVA) technique with 95 % confidence level by MINITAB V.16.1.It is important to know that for using ANOVA, three main hypotheses, normality, equality of variance and independence of residuals, must be examined.The results show that there is statistically significant difference between the obtained RPD values of the algorithms with the p values very close to zero for small and large sized problems.The results also shows that average of RPD for RGA is less than the other algorithms.Furthermore ach proposed algorithms, i.e.RGA and CVNS, are of better performance than the classical corresponding form, especially for large sized problems.

Conclusions and future researches
In this paper, a fuzzy multi-objective production distribution problem while unsatisfied demand is backordered was developed.The objectives are to minimize the total operational costs and to maximize the reliability of TS.The problem was mathematically formulated as a MILP.To examine the effectiveness of the model, different sized problems were designed and solved using two proposed algorithms called RGA and CVNS as well as the classical GA and VNS.RGA utilized several types of crossover operators and applied the better ones based on the historical performance.CVNS applied five different NSS.To control the algorithm, a learning procedure was designed.The algorithm recognizes which neighborhood performs well and tries to use those neighborhoods more than the others.The results indicate that RGA outperforms the other rival algorithms.Furthermore, RGA and CVNS outperform the corresponding classical forms i.e.GA and VNS for small and large sized problems.The proposed model can be useful for systems with delicate or sensitive products for which the type of TS is of high importance.Furthermore, it can be applicable for systems with very high shortage costs.One of the most important contributions of this paper is proposing two heuristic algorithms (i.e.RGA and CVNS) in order to solve the given model.Moreover, considering reliability for TS as a separate objective function is another contribution of this paper.
As future research, we can add supplier level to the current supply chain in order to provide raw materials required by the plants; this way, integration of TS with other operational aspects of the supply chain is extended.We can also consider a new objective function to minimize the delivery time to final customers.We can also assume uncertainty for other parameters of the model such as transportation time.

Fig. 2 Fig. 3
Fig. 2 The examples of both crossover operators.a The used one-point crossover operator.b The used multi-point crossover operator

Fig. 5
Fig. 5 The flowchart of the proposed CVNS

Fig. 6
Fig. 6 The convergence chart of the algorithms a GA, b RGA, c VNS, d CVNS

Fig. 7
Fig. 7 Average and interval RPD values of meta-heuristic algorithms for small and large sized problems.a Small sized instances, b Large sized instances To formulate the model, the following variables are required.The first five of these variables are continuous.They are concerned with the batch sizes of products, amounts of products transported by each TS from plants to DCs and from DCs to customers, amounts of inventories and amounts of delayed demand.The next two variables are of integer.They are concerned with the frequencies of vehicles traveling from plants to DCs and from DCs to customers.The last one is a binary variable to represent if each plant produces each product in each period.
2, 3, . . ., P) t Index of time periods (t = 1, 2, 3, . . ., T ) jkt The integer variable representing the frequencies of traveling of TS type m from DC j to customer k at time period t X pit The binary variable taking value 1 if plant i produces product p at time period t, 0 otherwise

Table 3
The optimal values of GA and RGA parameters for small and large sized instances