Ant Colony Optimization and Genetic Algorithm for Fuzzy Stochastic Production-Distribution Planning

: In this paper, a tactical Production-Distribution Planning (PDP) has been handled in a fuzzy and stochastic environment for supply chain systems (SCS) which has four echelons (suppliers, plants, warehouses, retailers) with multi-products, multi-transport paths, and multi-time periods. The mathematical model of fuzzy stochastic PDP is a NP-hard problem for large SCS because of the binary variables which determine the transportation paths between echelons of the SCS and cannot be solved by optimization packages. In this study, therefore, two new meta-heuristic algorithms have been developed for solving fuzzy stochastic PDP: Ant Colony Optimization (ACO) and Genetic Algorithm (GA). The proposed meta-heuristic algorithms are designed for route optimization in PDP and integrated with the GAMS optimization package in order to solve the remaining mathematical model which determines the other decisions in SCS, such as procurement decisions, production decisions, etc. The solution procedure in the literature has been extended by aggregating proposed meta-heuristic algorithms. The ACO and GA algorithms have been performed for test problems which are randomly generated. The results of the test problem showed that the both ACO and GA are capable to solve the NP-hard PDP for a big size SCS. However, GA produce better solutions than the ACO.


Introduction
In nowadays, as a result of globalization and technological developments, products have become more complex structure; the numbers and varieties of the components in products have increased, the number of suppliers has enlarged, and production systems have been transformed into a more complex system to be flexible. On the other side, distribution of products on time and at required quantity to the customer, who are located in different geographic areas in the world, became an important issue for the companies to improve the customer satisfaction level and profitability. Therefore, managers need the support of efficient approaches and methodologies to make right decisions in dynamic and complex market conditions. Production-Distribution Planning (PDP), which enables optimizing the supplying, manufacturing, and distributing process simultaneously, is one of the techniques to support the managers.
PDP aims to determine the varieties and quantities of raw materials that will be supplied and go on making decisions related to the production planning and distribution of final products to the customers [1]. It can be operated by two different approaches in supply chain systems (SCS); centralized and decentralized. In centralized approaches, a single member of the chain makes decision on PDP for the whole system. In contrary to the first approach, all the members of the SCS manage their own PDP in decentralized approach [2]. warehouses, retailers) with multi-products, multi-transport paths, and multi-time periods. The PDP handled in this study has been investigated by only Sakalli [1] in the literature according to our best knowledge. He developed a fuzzy stochastic 0-1 mixed integer mathematical model and converted it into a deterministic multi-objective mixed integer model by integrating fuzzy programming and chance-constrained programming techniques. Finally, he proposed a solution approach to solve deterministic multi-objective mix integer linear programming (MOMILP). The proposed solution approach can be successfully implemented into PDP for a small SCS. However, the mathematical model includes binary variables that select the transportation paths between the echelons of the SCS. Therefore, it cannot be solved optimally by using optimization packages when the numbers of the binary variables increase, which is a challenge.
The main contribution of this study is developing meta-heuristic algorithms to deal with this challenge. Two meta-heuristic algorithms have been proposed for big size PDP in a fuzzy and stochastic environment; ACO and GA. The proposed meta-heuristic algorithms have been designed for determining routing problem (selecting transportation paths) in PDP. The solution procedure mentioned by Sakalli [1] has been extended by integrating proposed meta-heuristic algorithms.
The paper is organized as follows: the problem definition and mathematical model are given in Section 2. Section 3 presents the proposed solution approach. The computational experiment is performed in Section 4. Finally, in Section 5, we make some concluding remarks.

Problem Formulation
The SCS for PDP includes four echelons. Plants manufacture products in regular and overtime by using raw materials supplied by different suppliers. The products are delivered to the customer by using warehouses and retailers. Customers only pick up their products from retailers. The SCS considered is drawn in Figure 1.  Sakalli [1] developed a 0-1 mixed integer programming model for PDP which includes uncertain parameters in both objective function and constraints. In order to solve uncertain model, it is required to transform the uncertain parameters into an equivalent deterministic one [51]. Therefore, he converted the fuzzy stochastic 0-1 mixed integer programming into a deterministic MOMILP by using fuzzy and chance-constraint approaches.
Fuzzy parameters, which are symbolized by "~", are modeled by using triangular fuzzy numbers. A triangular fuzzy number is denoted by the triplet (a p , a m , a o ) where a m represents the most possible value, a o represents the most optimistic value, and a p represents the most pessimistic value.
Fuzzy random parameters (transportation capacities), which are symbolized by "¯", are modeled by using the triangular fuzzy number as follows: Pr a , Pr b , and Pr c represent the probabilities of transportation capacity situations such as high, medium, and low capacities. Triangular fuzzy numbers (a p , a m , a o ), (b p , b m , b o ) and (c p , c m , c o ), define the amounts of each transportation.
There is one random fuzzy parameter CDP jrtc that represents the customer demand. The customer demand is a discrete variable which can occur in three situations such as high, medium, and low.
The probabilities of each situation are modeled as random fuzzy variables which are defined in Equation (2) for only high demand: ζ Pr(D) high = 0.8 with possibility 1.0 0.6 with possibility 0.8 On the other hand, demand quantities for each situation follow normal distributions with different fuzzy mean parameters and variances which is defined in Equation (3) for only high demand: where X θ high is a fuzzy variable. Therefore, demand quantity is a random fuzzy variable. Interested readers can be referred to Sakalli [1] for more detailed discussion on modeling uncertainties: The objective functions, Z 1 , Z 2 , and Z 3 define lower profit risk, possible profit, and higher profit possibility, respectively. Therefore, the decision maker aims to maximize Z 1 and Z 3 and minimize Z 2 . Equations (7) or (8) ensures that the total transported quantity from s for i at t will be less than, or equal to, on-hand quantity. Equations (9)-(11) are modeled for selecting the k from s to p and not to exceed its capacity where M is a big number. Equation (12) is a stock balance constraint for i in p at t. Equation (13) is a capacity constraint for inventory level of i in p. Equations (14)-(18) are capacity constraints for regular production, overtime production and product inventory levels in p, respectively. Equation (19) is an inventory balance constraint for j in p. Equations (20)- (22) are constructed to select the k from p to the w and not to exceed its capacity. Equation (23) is an inventory balance constraint in w. Equation (24) is an inventory capacity constraint for a w. Equations (25)- (27) are constructed to select k from w to r and not to exceed its capacity. Equation (28) is a balance constraint for inventory and backorder level in r. Equation (29) is an inventory capacity constraint for r. Equations (30) or (31) ensures to meet customer demand at a given probability level. Equation (32) is the definition of the decisions variables.
In order to solve PDP, the MOMILP is converted into an equivalent single-objective 0-1 mixed integer programming model by using Zimmermann's [52] fuzzy programming method. According to the fuzzy programming method, first, the positive ideal solutions (PIS) and negative ideal solutions (NIS) of the objective functions can be specified as Equation (33) [53]: Second, the linear membership function of each objective function is defined by Equations (34) and (35): Finally, the equivalent single-objective 0-1 mixed integer programming model is constructed as follows: maxλ Sakalli [1] proposed a procedure for the solution of PDP which is given as follows: Step 1: Construct the PDP model.
Step 2: Formulate the fuzzy and fuzzy random parameters as triangular fuzzy numbers and triangular fuzzy numbers with probability values, respectively. Model the demand quantities and situations as random fuzzy parameters.
Step 3: Transform the fully fuzzy objective function into three new crisp objective functions Step 4: Convert triangular fuzzy numbers into deterministic close intervals for a fixed α value. Use the lower bound of the close interval for constructing Equations (7), (14) and (16). Use the upper bound of the close interval for constructing Equations (8), (15) and (17).
Step 5: Convert the discrete fuzzy random parameters into deterministic close intervals. Use the lower bound of the close interval for constructing Equations (9), (20) and (25). Use the upper bound of the close interval for constructing Equations (10), (21) and (26).
Step 6: Convert the random fuzzy parameters into normally distributed random parameters with a deterministic close interval mean parameters Step 7: Fix an acceptable probability value (β). According to the β and chance-constraint approach, convert the demand constraint which is obtained at step 6, into deterministic linear constraints by considering lower and upper bounds of close interval separately.
Step 8: Find the maximum and minimum values for Z 1 , Z 2 , and Z 3 and construct the membership functions. Formulate a single-objective 0-1 mixed integer programming model (Z 4 ) by using membership functions. Step 9: Step 10: If the DM is not satisfied with the initial solution, return to Step 4 and update the α and β values and repeat the remaining steps until a satisfactory solution is found.
The proposed modeling and solution approaches can be solved global optimally by using GAMS's Cplex solver. However, if the numbers of echelons (s, p, w, r, k and t) increase, the number of binary variables in the model will increase exponentially. Therefore, optimization packages cannot solve the PDP global optimum, which is a reason to develop a meta-heuristic algorithm.

The Proposed Solution Algorithms
The 0-1 mixed integer PDP model is an NP-hard optimization problem for large SCS. In this study, the solution approaches of NP-hard PDP under fuzzy stochastic environments for a large SCS has been considered. Therefore, two meta-heuristic algorithms have been developed to solve PDP; the Ant Colony Algorithm (ACA) and Genetic Algorithm (GA).
The GA and ACO have been developed for route optimization process in PDP. The route optimization process aims to find optimal transportation paths between suppliers-plants (UKS), plants-warehouses (UKP), and warehouses-retailers (UKW). At each iteration of the GA and ACO, the obtained routes are sent to the GAMS as parameters and GAMS solves the remaining problem with an objective function value which is returned to GA/ACO as a fitness function value. The framework of the proposed solution approach is given in Figure 2. In order to integrate the proposed solution meta-heuristic algorithms, solution procedure, proposed by Sakalli [1], has been extended by reorganizing Steps 8 and 9, which is given in Section 2, as follows: Step 8: Find the maximum and minimum values for Z 1 , Z 2 , and Z 3 by using ACO or GA and construct the membership functions. Formulate a single-objective 0-1 mixed integer programming model (Z 4 ) by using membership functions.
Step 9: Solve Z 4 by using the ACO or GA.

Ant Colony Optimization (ACO)
ACO is an approximate optimization approach, proposed by Colorni, Dorigo, and Maniezzo [54], that was inspired by the behavior of ants. Ants use pheromones for sharing information in order to reach food by using the shortest path. ACO has been successfully implemented to several combinatorial optimization problems, such as Travelling Salesman [55], vehicle routing [56], and scheduling problems [57].
Chang and et al. [42] proposed a mathematical model for integrated PDP and a solution algorithm using ant colony optimization. Cheng and et al. [43] developed a 0-1 integer mathematical model for PDP and ant colony optimization to solve the production part. Calvete and et al. [44] proposed an ant colony optimization approach for decentralized PDP.
At the beginning of the ACO process, whole constants, variables (iteration limit, tolerance etc.), and weight matrix are defined. A weight matrix which will be updated according to the choices of ants and pheromone amounts is defined by the size of input variables (suppliers, plants, retailers, and warehouses).
Iteration limits and ant population sizes are defined in compliance with the objective function type. The optimization process proceeds till the iteration limit is reached. Each step of the iteration begins with a transportation route determination. Routes are defined by the weight matrix in each step of the iteration. At the beginning of the process, first route is defined randomly because of equality of the weight matrix. The existing route is sent to GAMS solver module and the success rate or with another saying objective function fitness value (z j ) is calculated and saved. Calculated objective value is compared to the best fitness value. If a new fitness value is equal or better than the tolerance level, the pheromone amount (ph) is calculated and each active relation of the weight matrix is updated and rewarded using the calculated ph value. If the fitness value is worse than the tolerance level, the weight matrix is updated by using the calculated ph value and the selection probability of the existing route is decreased. If GAMS finds an infeasible solution, the optimization process carries on with previous objective function value (z j−1 ). ACO process pseudo code for maximization problem is given in Table 1.

Genetic Algorithm (GA)
GAs are nature-inspired evolutionary search algorithms based on natural selection and genetics [58]. GAs give promising performance results for optimizing complex, large-scaled, single or multi-objective engineering problems and provides an optimal or a near-optimal solution [59]. The multi-objective optimization process is similar to a single-objective operation. However, in the single-objective process, there is only one global optimal solution; in the multi-objective process, there are more than one superior to the rest of the solutions in the search area [60]. Fundamentals of GA were developed by Holland [61]. Goldberg et al. presented the detailed research and implementation of GAs in various fields [62,63].
In a GA, each solution is called individual which is one of the global candidate solutions. The population consists of individual solutions and each individual carries ancestors' genetic information in their chromosomes.
The initial population is generally created with the randomly distributed dataset. The GA uses crossover and mutation operators to find the optimal solution. The crossover operator is used for information transference to new individuals and populations. Genetic diversity is held with the mutation operation.

Selection Operator
Selection operation is based on the idea of the "survival of the fittest". In general, the selection is the process of transferring the characteristics of the most powerful or individual having the closest fitness value to the next generation. There are many different selection methods in the literature such as roulette wheel, tournament, or truncation selection.

Crossover Operator
Global search called exploration is held with a crossover operator. Every individual consists of gene and chromosome sets and these genes and chromosomes carry over meaningful information coming from their ancestors. This operator chooses a crossing position or layout to diversify genetic pool. New individuals are produced by using this crossing combination. Breeding by using multi-point crossing is given in Figure 3.

Mutation Operator
The mutation operator is responsible for abstaining local optimal solutions and helps to prevent the recessive population. This operator randomly changes one or more genes values 1 to 0 or 0 to 1. The mutation constant is defined at the beginning of the optimization process.
The genetic algorithm-based optimization process begins with the transportation path (tp) relation base matrix creation. Each tp matrix involves a combination of three different relation matrices, which are UKS, UKP, and UKW. The initial population is involving four different randomly uniform distributed individuals. After the base population creation phase, all of the transportation paths (tp 1 to tp 4 ) are evaluated by using GAMS optimization package, respectively, and objective values are obtained. The roulette wheel selection method is used for elitism and provides transferring more powerful individuals' features to the next generations. Selected paths are used for generating new individuals by using single point cross-over operation. The crossing point is defined as the middle position of the path. After the crossing operation, each bit of the new individuals is controlled with the mutation constant. The mutation constant is defined as 1% before the optimization process at the beginning of the script. If randomly produced variable for each bit is greater than mutation value, it is changed to 1 to 0 or 0 to 1 randomly. After the mutation process, a new population is created and the optimization process continues while the iteration limit or deserved objective value is not reached. The optimization script is evaluated several times with 500-1000 iteration. The GA process pseudo-code for maximization problem is given in Table 2.

Computational Experiments
The computational experiments described in this section have been conducted in order to analyze the performance of the proposed solution procedures. According to our best knowledge, this is the first study that handles PDP for a large SCS which includes multiple-supplier, multiple-product, multiple-plant, multiple-warehouse, multiple-retailers, multiple-transport paths, and multiple time periods. Thus, any data for comparison does not exist. Therefore, five test problems have been considered to evaluate the proposed algorithms which are randomly generated.
On the other side, proposed ACO and GA have been compared with a Simulated Annealing (SA), which is a well-known and widely used algorithm in order to analyze and observe the superiorities (or weakness) of them. SA is a general form of optimization which refers to tempering materials such as alloys of metal, glass etc. by heating and cooling to reach perfect atomic structure [64]. SA is based on finding global optima using hill climbing method. Firstly, movement (position) is selected randomly from a pool of all possible movements and applied to the objective function. Then, if the current movement improves the solution it is accepted and the search area is moved to this position. Otherwise, the next movement is chosen from other ones, which are probably less than current [65]. The pseudo-code of SA is given in Table A5 in Appendix A.
In each test problem, there are six parameters: the number of suppliers (s), the number of plants (p), the number of warehouses (w), the number of retailers (r), the number of time periods (t), and the number of transportation paths (k) between echelons of the SCS. The characteristics of test problems are given in Table 3. There are five raw materials and one product in each test problem.  1  5  1  5  5  3  3  2  5  1  5  10  3  3  3  5  1  10  10  3  3  4  5  2  10  10  3  3  5  10  2  10  10  3  3 The proposed algorithms were implemented in MATLAB (The MathWorks Inc., New York, NY, USA), using the COIN-OR (Computational Infrastructure for Operations Research) solver of GAMS (GAMS Software GmbH, Frechen, Germany) optimization package. The computational experiments and validation were performed on a PC with Intel ®Core i7-7700K CPU at 4.20 GHz having 64 GB RAM.
In each test problem, α and β values were set as 0.4 and 0.95 respectively. Test problems were solved by COIN-OR and the optimal solutions for NIS and PIS of Z 1 , Z 2 , Z 3 , and Z 4 are given in rows which are entitled "Optimal" in Table 4. On the other hand, the proposed ACO and GA were repeated five times with 250 iterations for each test problem. The best solutions for NIS and PIS of Z 1 , Z 2 , Z 3 , and Z 4 obtained by the proposed approaches are reported in Table 4. The performances of the proposed algorithms are measured in terms of the solution gap (in %). The proposed algorithms achieved to obtain optimal solutions for Problem 1 and 2 which are small in size. When the size of the problem increases, such as in Problems 3, 4, and 5, the proposed approaches are capable of obtaining near-optimal solutions in an acceptable gap. It can be stated that the results of the ACO and GA approaches are close to each other while the GA finds relatively better solutions.
According to the gaps (%) given in Table 4, ACO and GA obtained lower values than SA in test Problems 3, 4, and 5. It is observed that the proposed ACO and GA approaches are more promising methods for fuzzy stochastic PDP than SA.
After the validation and verification of proposed approaches, a new test problem (Problem 6), which has a larger size than the other problems, was generated. Only ACO and GA have been performed to Problem 6 because of their robustness. Problem 6 includes 40 suppliers (s = 40) 12 raw materials (i = 12), three plants (p = 3), three products (j = 3), 12 warehouses (w = 12), 40 retailers (r = 40), and four transportation paths between each of echelons (supplier-plant, plant-warehouse, warehouse-retailer). Therefore, Problem 6 is more suitable for real-life problems. The parameters in the mathematical model which were randomly generated by using domains are given in Tables A1-A4 in Appendix A. Problem 6 could not be solved by using COIN-OR within the maximum limits of iteration and time. Z1 PIS , Z2 NIS , and Z3 NIS are minimization problems. The minimization process generally lasts longer than maximization process because of GAMS' solving procedures for the proposed methods and problems. Therefore, the iteration process is limited to 500 and 1000 for minimization and maximization problems, respectively. Comparative results of the proposed meta-heuristic algorithms are represented in Figure 4a-c Figures 5a-c and 6.   According to the negative and positive ideal solutions for Z 1 , Z 2 , and Z 3 , it can be observed that the GA shows better performance than ACO in four of six problems. However, these solutions are used as lower and upper limits for the membership function in Z 4 . The optimal solution of Z 4 cannot be obtained at the lower or upper bounds of membership functions because of conflicting objectives. Therefore, the comparative solutions of ACO and GA for Z 1 , Z 2 , and Z 3 are not significant indicators in determining the best meta-heuristic approach. The main problem that determines the best performance for proposed meta-heuristic algorithms is Z 4 .
According to the ACO, the positive ideal solution and negative ideal solution values of the objective functions are found as (15,742,124, 4,421,268), (45,856,319, 78,583,274), and (4,964,094, 15,906,761) for Z 1 , Z 2 , and Z 3 , respectively. When the equivalent single-objective 0-1 mixed integer programming model of the auxiliary MOMILP problem is solved, the total profit is obtained as a triangular possibility distribution with (54,920,893,66,474,301,78,332,275) and the overall degree of decision maker (DM) satisfaction with multiple goal values is achieved at 0.63. On the other hand, GA obtained the positive ideal solution and negative ideal solution values of the objective functions as (15,757,821,4,694,713), (45,781,117,79,004,432), and (4,852,275, 15,781,522) for Z 1 , Z 2 , and Z 3 , respectively. When the equivalent single-objective 0-1 mixed integer programming model is solved, the total profit is obtained as a triangular possibility distribution with (56,243,986,68,505,864,80,833,744) and the overall degree of DM satisfaction with multiple goal values is achieved at 0.684., Total profits in triangular fuzzy numbers can be defuzzified by using the expected value (EV) of fuzzy numbers which is given in Equation (37): The expected profit values of the ACO and GA have been calculated as 66,550,442 and 68,522,364 Turkish Lira, respectively, by using Equation (37). According to the total profit and DM satisfaction level, it can be concluded that GA found better solutions than ACO.
The solutions of the proposed algorithms are summarized in Tables 5-7 because of large sizes of indices. Although ACO produces more product 1 and 2 (j1 and j2), its production quantity for product 3 is significantly smaller than the GA. Therefore, GA produced and sold more products than the ACO in total. The available capacities of the plants are not adequate to meet all demand. This situation caused to backorders. Therefore, there are no inventories in warehouses and retailers. Total backorder level of GA is smaller than ACO which enables reducing the cost in GA.
According to the total purchased raw material quantities by plants, it can be absorbed that the numbers of raw materials supplied GA and ACO are close to each other. However, while GA supplied raw materials from 13 different suppliers, ACO used 23 suppliers. Cooperating with more suppliers means extra fixed transportation.  Table 7. Summary of results.

Conclusions
In this paper, a fuzzy stochastic PDP has been considered for an SCS which includes four echelons. The PDP handled in this study was considered by Sakalli [1] and modeled as a 0-1 mixed integer model. He developed a solution procedure and successfully implemented it for small SCS. However, it is not possible to solve the mathematical model global optimally by using optimization packages for large SCS because of the binary variables which are designed for selecting routes between echelons in SCS. Therefore, it is required to develop meta-heuristic algorithms to solve NP-hard PDP.
In this study, we have developed two meta-heuristic algorithms, GA and ACO, in order to solve Sakalli's model for large size SCS. In the solution process of PDP, meta-heuristic algorithms (GA/ACO) are performed for route optimization. By determining routes, 0-1 mixed integer model is converted into a deterministic linear programming model and can be solved optimally by the COIN-OR solver of GAMS.
The proposed solution approaches have been performed for randomly generated test problems. The results of the test problems have showed that both proposed meta-heuristic algorithms are capable of solving the problem where the GA obtain better solutions than the ACO.
In real-life applications, some parameters of the PDP can be affected by the conditions of dynamic market structure and it is not possible to define them precisely in the decision making process at the tactical or strategic level. This challenge can be overcome by including the manager's expertise and judgments into the modeling process by using fuzzy, random fuzzy, and fuzzy random parameters. Therefore, decision or policy makers can manage the large production and distribution systems more effectively by using the proposed approach. The proposed approach can be used in several industries, such as automotive, electronic devices, and textiles, by implementing the proposed extended procedure in ten steps. The proposed approach is implemented iteratively for real-life problems. In each solution, the proposed approach produces several performance values that show capacity usage proportions in the production-distribution system. In this way, decision makers can detect the bottlenecks in the system and take preventive measures for balancing the capacities of supply, transport, manufacture, and storage. On the other hand, the proposed approach can be used to manage customer demand. The proposed approach enables the analysis of the results of changes in demand quantity on profitability by using random fuzzy parameters. Therefore, managers can specify their marketing and pricing strategies. Consequently, decision makers can evaluate the effects of different predictions and judgments by performing the proposed approach several times.
There are two main limitations of the proposed approach. The first one is determining the minimum acceptable possibility degree (α) which represents the decision maker satisfaction level and takes values between 0 and 1. A very low possibility degree may cause obtaining an optimistic and unrealistic solution. On the contrary, a very high possibility degree may cause obtaining a conservative solution which can be feasible, but not robust. Therefore, the proposed approach has to be performed with different α values in real-life applications, which means longer processing time. The second one is starting with the GA/ACO with a better feasible solution. According to the nature of proposed approaches, the starting solutions have just been selected randomly. This situation directly affects the processing time which is required to approximate the optimal solution. If GA/ACO start from a better feasible solution, they reach the best solution which is close to the optimal value in a shorter processing time.
With respect to future research, it is possible to integrate other heuristic algorithms, such as a Flower Pollination Algorithm (FPA), Krill Herd Algorithm (KH), and Particle Swarm Optimization etc., into the GA/ACO in order to obtain better starting solutions. On the other hand, PDP can be modeled as bi-objective functions by considering the maximization of the customer satisfaction level or minimization of the total transportation time.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The MOMILP is given at below. holding cost of j in p at t RHC j required storing capacity for j TCPW pwkt transportation capacity of k from p to w at t FTCP pwkt fixed cost of using k from p to w at t VTCP jpwkt variable cost of transporting j from p to w by using k at t RTC j required transporting capacity of j TCWR wrkt transportation capacity of k from w to r at t FTCW wrkt fixed cost of using k from w to r at t VTCW jwrkt variable cost of transporting j from w to r by using k at t W HC jwt holding cost of j in w at t W IC w products storage capacity in w HCR jrt holding cost of j in r at t RIC r storage capacity in r BCR jrt backorder cost of j in r at t CDP jrt demand of c for j from r at t POP jt price of j at period t Decision Variables TRQ ispkt transported quantity of i from s to p by using k at t SRP ipt stored quantity of i in p at the end of t RPQ jpt quantity of j produced at regular time in p at t OPQ jpt quantity of j produced at over time in part TPQP jpwkt transported quantity of j from p to w by using k at t SLP jpt stored quantity of j in p at the end of t SLW jwt stored quantity of j in w at the end of t TPQW jwrkt transported quantity of j from w to r by using k at t SLR jrt stored quantity of j in r at the end of t BLR jrtc backorder quantity of j in r for c at the end of t SPQ jrtc sold quantity of j from r to c at t UKS spkt 1 i f k is used between s and p at t 0 otherwise UKP pwkt 1 i f k is used between p and w at t 0 otherwise UKW wrkt 1 i f k is used between w and r at t 0 otherwise Appendix A   Table A5. Pseudo-code for SA.