Improved Meta-Heuristic for Multi-Objective Integrated Procurement, Production and Distribution Problem of Supply Chain Network Under Fuzziness Uncertainties

: In this study, we devoted a design under uncertainty of a supply chain network composed of multiple suppliers, multiple plants, multiple distributors and multiple customers. The proposed model is a bi- objective mixed integer linear programming which consider several constraints and optimize the total costs including the procurement, production, storage and distribution costs as well as maximize on-time deliveries. In order to model a similar problem to the real planning problems, the objective function coefficients (e.g., procurement cost, production cost and inventory holding) and other parameters (e.g., demand and production capacity), are all supposed triangular fuzzy numbers. Besides, a hybrid credibility-based model is constructed for the problem, i.e., expected value and chance constrained models. Moreover, the construction of the crisp equivalent model is based on different properties of the credibility theory. The resulted crisp equivalent model is a bi-objective MILP. To transform the crisp BOMILP to a single objective MILP model, we propose the LP-metric method. To solve the planned model for largest instances, we tend to propose a hybrid algorithm based on a recent metaheuristic called Hunger Game Search (HGS). This algorithm employs the so-called chaotic search to prevent premature convergence and entrapment in local optima. Finally, the efficiency and sensitivity of the proposed model are justified through the numerical results.


Introduction
Given the current trends in globalization Douaioui et al. (2018), many companies are increasingly sourcing, producing and marketing their products around the world. For many of them, success depends on a personalized production, in a short time and at a good price. Recently, supply chain management has been proven as the best solution to meet these contradictory constraints. As a result, many researchers focus on planning and scheduling of supply chain in order to maximize the overall performance. This integration not only reduces the number of steps in the process but also tends to eliminate the inherent barriers between the various functions to see overall optimization Douaioui et al. (2021). Our research is in line with this trend, intending to study the problem of integration of functions such as procurement, production and distribution of the supply chain. It should also be noted that the challenges associated with the variability of demand are accentuated by the length of the manufacturing period and the transport conditions. All these factors make planning and scheduling of the supply chain a challenged task. The problems cited above are often dealt with under specific conditions. As reported by Klibi et al. (2010), any supply chain planning that relies on deterministic conditions risks losing its durability. They also mention that in some cases, it is not enough for the company to consider usual parameters such as demand, prices, or other parameters such as random variables, but undesirable events such as terrorist attacks and natural disasters. Consequently, it's necessary to develop a adequate strategy that integrates uncertainty into supply chain network.

Supply Chain Network
This study highlights the problem of integrating functions of supply chain namely procurement, production and distribution and which has several geographically dispersed suppliers, production sites, distribution centers and customers. Each supplier has a limited procurement capacity. Each production has a limited production time capacity which depends essentially on the installed processing units and the speed of each unit. The utilization of a production unit generates production and setup costs that must be integrated into the model. In general, the optimization procedure must consider the other costs such as holding inventory, transportation and backorder. The model considers also the on-time delivery to minimize too late the product of each customer. Figure 1 and 2 represents an illustration of the supply chain network. Regarding the significant increase in supply chain complexities, several researchers developed and solved models by the metaheuristic approach. Among them, Samadi et al. (2018) develop a multi-objective closed-loop supply chain, which has considered the total cost, environmental impacts of the network and the social objective. Abdi et al. (2021) address a new closed-loop supply chain model using the two-stage stochastic programming formulation. The authors solve the model by Whale Optimization Algorithm (WOA) and Particle Swarm Optimization (PSO), the used meta heuristic performance outperformed the Genetic Algorithm (GA) and Simulated Annealing (SA). Bottani et al. (2019) propose a new resilient food supply chain design problem, the authors propose a bi-objective mixed-integer programming formulation. The paper aims to minimize the lead time and maximize the total profit. The authors use Ant Colony Optimization (ACO) to solve the problem in a large-scale network. Sarkar et al. (2019) address a multi-attribute closed-loop supply chain model to maximize total profit and minimize the emissions of CO2. The paper uses three metaheuristics; Genetic algorithm, Interior point optimization algorithm and Particle swarm optimization to solve the model in large size instances. Khalifehzadeh et al. (2019) present a multi-echelon production-distribution network model in the multi-time horizon, the model aims to minimize the total delivery lead time and maximizes the total profit. The proposed solution shows that the performance is improved by applying the Firefly Algorithm and Ranking Genetic Algorithm, respectively. Regarding complexity and the multiplication of interevent in supply chain network, several researchers like this study contribute the new meta-heuristic to improve the overall performance.

Supply Chain Management Under Uncertainty
The theory of probability proved to be the old theory to deal with uncertain parameters. However, in practice, several disadvantages are related to the use of the probability theory (1) need for a sufficient and reliable history and (2) problem of modeling of subjective parameters Pishvaee and Torabi (2010), Liu and Kao (2004). The fuzzy set theory (FST) is an alternative to treat uncertainty by giving local preferences into account in optimization problems. In the fuzzy set theory, there are three forms of measures for a fuzzy event namely the possibility, the necessity and the credibility measure. Possibility theory is considered as a counterpart of probability theory Dubois et al. (2004) and is widely used in literature to treat fuzzy variables. Further, it is conflicting with the law of excluded middle and the law of contradiction. To over this challenge, (Liu and Liu, 2002) proposed the credibility measure. This measure is a more reasonable fuzzy inequality indicator than the possibility and necessity because it compensates for their disadvantages. For example, an event with a maximum possibility of 1 might not happen while an event with maximum credibility of 1 occurs. Furthermore, a fuzzy event with maximum possibility 1 sometimes carries no information while a fuzzy event with maximum credibility 1 means that the event will happen at the greatest chance Huang (2006). In the optimization process for management and planning, it is usually assumed that the credibility level should be super than 0.5 in response to avoiding un-satisfactions and overtaking the risks Pishvaee et al. (2012).X supposed be a fuzzy variable with membership function m(x) and t a real number. Based on (Liu and Liu, 2002) the credibility measure is characterized as Eq. 1: Appropriately, the credibility measure can0 be expressed as an average of the possibility (Pos) and necessity (Nec) measures. In addition, the expected value of  is characterized as Eq. 3 Liu and Liu (2002): Now, suppose that  is a triangular fuzzy number expressed by three points as  = (a, b, c Based on (4) and (5), it can be demonstrated Zhu and Zhang (2009) hat if  is a triangular fuzzy number and  > 0.5 correspondingly: Equation 6 and 7 can be performed and a lot of conveniently when compared to a-critical values proposed by Liu and Kao (2004), to convert fuzzy chance constraints to the equivalent crisp model. In credibility programming, there are three forms of credibility approaches: The chance constrained programming, the expected value Liu and Liu (2002) and the dependent chance constrained programming Liu (1999). The chance constrained programming model uses the expected value operator for each imprecise coefficient. It can be applied easily without increasing the complexity of the original model compared to the other two methods, but at the same time, it has no control over the level of confidence of the fuzzy chance constraint. The second model can control the level of satisfaction of the fuzzy chance constraint by using the concept of b-level. Although increases the complexity of the model since it adds a new constraint for each objective function of the main model. In this study, we combined the expected value operator and the chance constraint to deal with the uncertain parameters of the objective functions and constraints. The reason behind this consolidation is that the expected value operator does not need any additional data for objective function or the ideal solution and is also benefited from advantages of the chance-constrained programming the approach to control the level of satisfaction of the fuzzy chance constraint by using the concept of α-level Pishvaee et al. (2012).

Mathematical Formulation
According to the above-mentioned descriptions and justifications, the proposed credibility based fuzzy mathematical programming with expected value operators for integrated production-distribution in the supply chain network. The basic considerations of the problem under consideration are summarized, as follows: The planning horizon is known and divided into a set of period t ∈T A set of suppliers s ∈S, production sites p ∈ P, a set of DCs d ∈ D and a set of customers c ∈ C A set of products k ∈ K.
The key decision variables are:  The set-up of products to units in each period t  The quantity produced of each product on each processing unit at each period t  The total quantity dispatched between production sites-DCs and between DCs-Customers respectively  The inventory levels of each product at sites and distribution center DCs at the end of each period  The backorder quantity of each product p provided by the customer at each period

Indices and Sets
 s: Index used for a supplier s = 1,…, S  r: Index used for a raw material r = 1,…, R  p: Index used for a plant p = 1,…, P  k: Index used for a finished product k = 1,…, K  d: Index used for a Distribution Center (DC) d = 1,…, D  c: Index used for a customer c = 1,…, C  t: Index used for a planning period t = 1,…,T

Decision Variables
 Xk,p,t: 1-if product k is produced by plant p in period t  IRPr,p,t: Inventory level of product r in plant p at the end of period t  IKPk,p,t: Inventory level of product k in plant p at the end of period t  IKDk,d,t: Inventory level of product k in DC d at the end of period t.  QRSPr,s,p,t: Quantity of raw material r dispatched from supplier s to plant p in period t  QKPk,p,t: Quantity of finished product k produced in plant p in period t  QKPDk,p,d,t: Quantity of finished product k dispatched from plant p to DC d in period t  QKDCr,d,c,t: Quantity of finished product k dispatched from DC d to customer c in period t  BQKCk,c,t: Quantity of backorder for product k incurred by customer c in period t

Objective Functions
We consider two important objectives in our supply chain network problem: The first objective (TC) is the principal aim of the major companies. The objective is to minimize the cost in different operation of supply chain and to maximize the exploitation servals resources. It consists mainly of the total costs such as: Procurement costs, Inventory costs, transportation costs and Backorder costs.
Pr, , , * P r, , * P The second Objective function (OTD) is considered to maximize the on-time deliveries to ensure a continuous flow of production and minimize the plant stoppages due to unavailability of raw materials items. In this study, we have assumed that the demands, supplier capacity and all units' costs (i.e., production, transportation, storage and backorder are all fuzzy variables. So, the total costs and the on-time deliveries of raw materials are also fuzzy variables. In the process of search an optimal solution, it's not sufficient to maximize on-time deliveries and minimize the total costs, but it's primordial to specify the fuzzy variable. Therefore, we develop a new model within the framework of credibility theory given as follows:

Constraints
Constraint 10 granted that the total quantity of raw materials that can be transported from each available supplier to any manufacturer cannot surpass the fuzzy supplier capacity: r, , , P , , ; , Constraints (11)(12) give the inventory balance equations of raw materials and finished products in plants respectively. r, , 1 r, , , , , , , Constraint Eq. 13 ensure the inventory balance equations of finished products in DCs:   , ,  , , 1  , , , , , , ; , , Constraint Eq. 14 guarantees that the fuzzy total required time to produce the products cannot exceed the fuzzy total available time: , Constraint Eq. 15 shows the quantity of the raw material shipped to plant p from suppliers plus the inventory level of raw material in plant p in period t is limited to the capacity of plant: , Constraint Eq. 16 assure that the volume of production must be less than or equal to the total storage capacity of the plant p in time period t: , , , Constraint Eq. 17 granted that the quantity of each product shipped to each distribution center from any plant plus the inventory level of product k in the distribution center d in period t is limited by the capacity storage of the distribution center: Constraint Eq. 18 represents that the backlog level. In the last period t equal to backlog level of the past period plus the expected value of demand minus the total received from DCs: , , Constraint 21 guarantees that the backorder quantities of product k in the distribution centers (DC) d in time period t is limited to a percentage of the demand of that period: Constraint 22 guarantees that the on-time delivery of finished products from DCs to customer c: , , Constraint 23 defines the status of the decision variables:

Transforming a Multi-Objective Problem into a Single Objective
Generally, there is a difference between the singleobjective meta-heuristics and multi-objective metaheuristics: The multi-objective meta-heuristics provides a compromise between contradictory objective functions that result in a set of optimal non-dominated solutions called Pareto-optimal solutions, while the single-objective metaheuristics upon obtaining a single optimal solution are proved. Any objective can be modified only at the cost of ruining another objective/s according to each Pareto optimal solution. So, for all that, the decision-makers try to articulate their posterior preference or prior over the objectives based on multi-objective techniques. Thus, the scalarization is used to transform a multi-objective problems to single-objective problems. In this study to solve the mathematical model, the Compromise Programming (CP) method is considered. The aim of this method is to search a solution that comes as close as possible to optimal (ideal) values of each objective function. The closeness is defined by L-p distance metric as follows: In which f1.f2… fk are different and conflicting objective functions and * * * 12 ... k f f f are the ideal value for the i th objective. The largest deviation according to the increasing the is dedicated. Moreover, the importance of this method can be justified by the detection of the proper value for r; the value that fits the attitude of the decision-maker). According to the priority of expressed decision-makers, this method can generate only a single Pareto solution. In this study, r is chosen 1. Since the model proposed is strongly NP-hard a solving methodology, based on GA is adopted to resolve the model in the next section.

Fuzzy Simulation
Since total measure costs TC (Q; ξ) is a fuzzy variable, it will be difficult to the expected value by using analytical methods. As a kind of Monte Carlo methods, fuzzy simulation (Selim and Ozkarahan, 2006) provides an effective approximation. The main steps of fuzzy simulation are given according to (Liu and Liu, 2009).

Brief Description of Hunger Games Search (HGS) Optimization
The Hunger Games Search (HGS) is a new population-based model to deal with optimization problems based on social animals' characteristics in searching for food. Precisely, the algorithm hunger games search turns as follows: In each iteration, the Algorithm (HGS) searches near the optimal location, where the weights values, or the hunger values, represent the influence of hunger on the animal's individual activities. The adaptive and time-varying mechanisms allow to hunger games search algorithm the ability to deal with multi-modality and local optimal problems more efficaciously.

Steps of Fuzzy simulation
Step 1. Set e = 0.
Step 7. Repeat the fourth to sixth steps for N times.

Approach Food
To express its approaching behaviors in mathematical formulas, the following formulas are proposed to imitate the contraction mode: where, R is in the range of [-a, a]; r1 and r2 respectively represent random numbers, which are in the range of [0, 1]; randn (1) is a random number satisfying normal distribution; t indicates that the current iterations; 1 W and 2 W represent the weights of hunger; b X represents the location information of a random individual in all the optimal individuals;

 
Xt represents each individual's location; and the value of l has been discussed in the parameter setting experiment.
The formula of E is as follows: where, i  1,2,…, n, F(i) represents the fitness value of each individual; and BF is the best fitness obtained in the current iteration process. Sech is a hyperbolic function The formula of R is as follows: where, rand is a random number in the range of [0, 1]; and Max_iter stands for the largest number of iterations.

Hunger Role
The starvation characteristics of individuals in search are simulated mathematically.
where, AllFitness (i) preserves the fitness of each individual in the current iteration.
The formula for H can be seen as follows: Represents the random number in the range of [0,1] F(i): Represents the fitness value of each individual BF: Represents the best fitness obtained in the current iteration process WF: Represents the stands for the worst fitness obtained in the current iteration process UB and LB: Indicates the upper and lower bounds of the search space, respectively H: Represents the hunger sensation, is limited to a lower bound, LH

Brief Description of Chaotic Search
The chaos process is sensitive to the initial condition Schuster (2005). The nature chaos process is deterministic and quasi-random. So, the chaos process is apparently random and unpredictable. Mathematically, the chaos process is a randomness of a simple deterministic dynamical system. Then, the chaotic map can be defined as discrete-time of dynamical system running in a chaotic condition Alatas (2010): where, {xk :1, k = 0,1,2,…,N} : Represents the chaotic sequence, which can be utilized as spread-spectrum sequence as random number sequence. The chaotic sequences have been confirmed to be simple and fast produce; it is needless to store long sequences Heidari-Bateni and McGillem (1994). Only a few functions (chaotic maps) and parameters (initial conditions) are required even for very long sequences Alatas (2010). In this study, r chaotic variables are generated by the logistic mapping technique represented in Eq. 34; in this study we use the logistic map for four reasons:  Firstly, the new position is found out by the logistic map and the jump over local optimum will not be repeated. So, the speed of convergence can be accelerated  Secondly, the logistic map is simple and so very easy to implement  Thirdly, the logistic map is easy to fix in every part of the hunger games search algorithm  Fourthly, the complexity of the hunger games search algorithm will increase since chaotic map, provides every dimension, using another chaotic map: Let, j = 1,2,…, N-1 and hence, other N-1 solutions are produced by the same method.

Solution Representation
When using hunger games search algorithm an importing task is how to convert the solution of a problem into agent code. It is necessary to analyze all application problems by combining the expression of the solution and the evolution operator that fits the problem. In the literature many encoding schemas are proposed Yu et al. (2006). In this study, the structure of problem's agent includes binary and integer variables which consist of six parts. The first part of agent indicates the procurement quantity between supplier and plant in period t. The second part is a binary which indicates if the production of product take place in plant or not in period t. The third part indicates the produced quantity in period t, while four and five part indicates the transported quantity between plants-DCs and DCs-Customers respectively. Finally, the six parts provides the incurred back-order quantity from product from customer in period t.

Fitness Function
This study deals with the fitness function based on sequence. Np agents of the population are arranged in descending order according to the values of the target function, which means the agents are arranged from good to bad. Suppose the parameter a (0;1) is given then the fitness function Liu and Liu (2009) based on sequence is: where i = one demonstrates that the agent is the best solution whilei 1 /4 pop size; reveal the worst. As a result the best agent does not exist in the last generation., the best agent marked V0 has to be reserved at the beginning calculation of iteration. If there is a better agent in the next generation then we use it to replace the former V0. At the end of evolution, this agent can be considered as the solution of the model.

Execution of Test Termination
The stop criterion is determined according to the parameter max gen. If iter < Maxgen, continue with Step 3 to Step 6 in the hybrid intelligence algorithm process. Ifiter > Max gen, stop the algorithm. The maximum fitness value and its corresponding individual of the population in the total iterative process are recorded.

Computational Results
In this study, the problem of integrating the procurement, production and distribution functions of a four-echelon supply chain network is presented, Table 1 shows the related random distribution for each parameter. In this study, we use the triangular fuzzy numbers to model imprecise parameters rather than others fuzzy distributions such as trapezoidal due to their simplicity and computational efficiency.
The fuzzy numbers ξ = (ξ p , ξ m ; ξ o ); where, ξ p is the most pessimistic, ξm is the most likely and ξ o is most optimistic values. These values must be estimated for each fuzzy parameter.  In doing so, the method proposed by Lai and Hwang (1993) is used. First, the most likely xim value for each imprecise parameter is specified randomly corresponding to the uniform distribution, the used values are illustrated in in Table 1. Then, the most pessimistic ξ p and the most optimistic ξ o values of a fuzzy number ξ are obtained as ξ p = (1 -r1) ξ m , ξ o = (1 + r2) ξ m where (r1, r2) are two numbers randomly generated according to the uniform distribution [0.1,0.3]. The predetermined level is set to 0.8. The sizes of the designed test problems are given in Table 2.
The equivalent auxiliary crisp model is coded in GAMS 22.5/CPLEX 12.2 software and all experiments are computed using a Core i5 2.10 GHz computer with 16 GB RAM.

Sensitivity and Robustness Analysis
It follows from Fig. 3 that the storage costs including (storage in plants and DCs) are equal to zero due to the fact of the favor of the distribution (distribution costs represent the highest value in the graph) rather than storage to satisfy the customer demand which verified by the very small backlog costs rather than production or distribution costs. In addition, we carried out the sensitivity analysis in order to investigate the impact of some input parameters affecting the both objective functions such as the capacity of supplier, the maximum of capacity of production. Table 3 represents the impact of on-time delivery provided by suppliers as well as the supplier capacity on values of the second objective. It follows from Fig. 4 that the quantities of raw materials delivered on time increase as both parameters increase. A decision is to consider the highest weight to these criteria when selecting suppliers.
In all benchmarks, we use Standard Deviation (STD) and the Average (AVE) to compare all metaheuristics. The using of these two measures demonstrate the ability of algorithms to avoid local minima. The Table 4 demonstrate that the average objective values found by HGS and CSHGS over perform the ones obtained by the GA and PSO in different instances. The convergence curve of the different metaheuristic in optimizing problem instance 1, 3 and 6 are given in Fig. 5-9. The instances are divided in two types of problems as small size and large size. In the small size, to ensure the integrity and accuracy of the model, their results are obtained by programming mathematical model in GAMS software. Table 3 demonstrates an objective function value for each problem with various indicators and parameters in small size. In the large size, four test problems are examined that show the results of the proposed method compared with the well-known methods such GA, PSO and the standard HGS. These comparisons are calculated based on the deviation of the current value and best value found so far. Tables 4-5 shows computational results which involve combined objective function in the small and large sizes. Moreover, the convergence performance is given of all algorithms of some instances.      Instance 1 2  3  2  3  2  2  10  Instance 2 4  5  4  5  4  3  10  Instance 3 6  7  6  7  6  6  10  Instance 4 8  9  8  9  8  8  10  Instance 5 10  11  10  11  10  10  10  Instance 6 12  13  12  13  12  12 10       --k2  p3  1063,525  -----p4  1219,026  -----p7  ----4372,894  16510,461  p2  2315,674  2315,674  ----k3  p4  -----13446,979  p5  --1286,509  ---p6  -6614,733  ----k4 p2 ---6523,538 5204,306 5204,306 As it is represented, the GA and PSO get stuck into local optimal solutions after a while iteration, respectively while the CSHGS explore the space search in early iteration and can reach convergence in later iterations with the best identified result 0.208. Consequently, the solution of the CSHGS outperforms the other algorithms in search quality and uniform over iterations. The optimal solutions are presented in Table 6-8.

Conclusion
This study devoted the integration problem of procurement, production and distribution in supply chain network under uncertainty within the credibility framework. The model proposed in this study is a biobjective mixed-integer linear programming. This model considers technological constraints issuing a lot of supply chain network in world that aims to minimize the total costs such as production, storage and distribution as well as maximize the on-time deliveries to ensure a continuous flow of production and minimize line stoppages in plants due to unavailability of raw materials items.. To model the problem which is similar to the real planning problems, the objective function coefficients (e.g., procurement cost, production cost and inventory holding) and other parameters (e.g., demand and production capacity). The proposed model combines the expected value model and chance-constrained model. To transform the bi-objective model to single objective model the Lp-metric is used. To solve the proposed model for largest instances, we address a hybrid algorithm based on a recent metaheuristic called Hunger Game Search (HGS). This algorithm employs the so-called chaotic search to prevent premature convergence and entrapment in local optima. The numerical results showing the efficiency of the proposed solution. The proposed model considered the fuzziness as the source of uncertainty. However, future studies should aim to address hybrid uncertainties, such as an encounter with fuzziness and roughness. For instance, it is broadly accepted the representation of demands a triangular fuzzy number (a, b, c) variable from the view point of the fuzzy theory, but the values of a, b and c may emerge within complete or uncertain information. In a sense, the values of a, b and c are of rough characteristics. Thus, the decision-makers have to manage the fuzzy number of rough parameters by adapting the clients demand in order to be fitly defined as a fuzzy rough variable Liu and Liu (2002).