Considering discrete delivery ordering and shipment consolidation in a three-echelon supply chain with failure

: The aim of this paper is to extend research results published in two recent papers which had investigated single-vendor single-buyer inventory model with lost sales and discrete delivery orders and a model in which a single buyer orders a product from multiple vendors who deliver their products to the buyer as joint shipments. We extend the two models studied in the literature by considering a multi-product three-echelon supply chain including multiple suppliers, single manufacturer and single distributor with discrete delivery orders for raw materials and products. We assume that each product needs a specific raw material provided by the suppliers and the raw materials are delivered in batches to the manufacturer, which produces n products through n types of raw materials, and sends them to the distributor. All the products are manufactured utilising one machine and delivered in batches to the distributor. We also assume that a portion of the batches sent to the distributor by the manufacturer are defective. We have developed two heuristics based on genetic algorithm (GA) and biogeography-based optimisation (BBO) algorithm and our numerical results indicate that the GA-based heuristic slightly outperforms that of the BBO algorithm.


Introduction
Inventory overhead is one of the costly elements in many organisations, but because its effective deployment contributes to business success, organisations invest in inventory management of raw materials, work in process, finished goods and etc. With globalisation trends and increase in competition, inventory management has become a key factor to remain in today's competitive markets since customers expect to receive their commodities quickly. Despite extensive researches in this field, there is a significant gap between real-world problems and existing theoretical researches.
Appropriate inventory management is one of the challenging activities of every production system since excessive inventory accumulation results in unnecessary holding costs, e.g., perishability, obsolescence, capital recession, insurance and etc. Therefore, determining optimal inventory level for accurately estimated production requirements is quite essential. Such approach can trigger a proper balance between customer demand and production capacity, reduce ordering, production and inventory costs and increases customers' satisfaction.
One of the key questions about firms involved with outside suppliers is to determine the order quantity and the reorder point. One of the first, simplest and widely used studies in inventory control problems is economic order quantity (EOQ) model proposed by Ford Whitman Harris (Cárdenas-Barrón et al., 2014a). He was the first researcher who investigated optimal lot sizing in production systems, by publishing a paper titled: "How many parts to make at once?" in February 1913. Harris assumed that orders are received at once; however in production systems this assumption is not practical. As a result, a modified version of this model was proposed by Taft (1918) by considering production rate and proposing basic economic production quantity (EPQ) model for manufacturing units.
Although, the EOQ and EPQ models are widely employed in the real-world problems (Tersine, 1994;Silver et al., 1998), their inherent assumptions for the addressed basic models restrain their actual implementation. Many researchers have tried to extend the addressed basic models in the past few decades.
In this paper, a multi-product three-echelon supply chain including multiple suppliers, single manufacturer and single distributor with discrete ordering system is considered. Each product needs a specific raw material provided by the suppliers. The raw materials are delivered in batches to the manufacturer. All the products are manufactured utilising a machine and delivered in batches to the distributor. Some practical examples for this kind of problem can be found in the manufacturing of concentrated juice and fruit packaging industries. Each kind of fruit can be purchased from an outside supplier (farmer) or own farm in batches (truck loads). However, considering seasonality of fresh fruit production, items can also be purchased from third parties who possess both farms and cold storages. The reason behind purchasing each fruit from an independent supplier is related to the dependency of fruit's quality to its nurturing climate.
The structure of the rest of the paper is as follows: literature review and mathematical model are presented in Sections 2 and 3, respectively. Section 4 describes solution procedure. Section 5 presents computational results. Conclusions and future research recommendations are given in Section 6.

Literature review
Batch production and shipment is one of the main concerns in production systems and inventory control models. This concern stems from real world cases where the system has items produced in batches with a known size, for instance production of bottle caps by cap compression molding machines and so forth. One of the first studies in this area has been performed by Pasandideh and Niaki (2008). They studied multi-product EPQ model considering discrete order delivery with storage constraint. All products are produced by a single supplier and given to a company in pallets. The company as the customer makes all the decisions on pallet sizes and the number of shipments. Furthermore, a minimal and maximal number of shipments are considered. In another study,  optimised discrete deliveries in EPQ inventory model, i.e., the contractor shipped orders in form of pallets, for a problem with one product where shortage and delays were not allowed.  extended Pasandideh and Niaki (2008) model considering unsatisfied demand to be backlogged. They proposed a parameter tuned genetic algorithm (GA) to solve the problem. Widyadanaa and Wee (2009) extended Pasandideh and Niaki (2008) as a mixed integer non-linear programming (MINLP) model with budget constraint under Just in Time (JIT) mode of operation. They developed a branch and bound solution procedure using Lagrangean technique in order to reach global optimum. Recently, Cárdenas- Barrón et al. (2014b) proposed two heuristic algorithms in order to solve the addressed problem. The solution time of the given heuristics increase linearly while that of Widyadanaa and Wee (2009) increases exponentially with the problem size.
Another assumption considered in production planning and inventory control systems is producing all products with same machine. Some of the first studies that studied the aforementioned assumption are Bomberger (1966), Madigan (1968), Baker (1970) and Elmaghraby (1978). Due to the complexity of single machine problems, no exact solution procedure for them could be provided. However, Hanssmann (1962) proposed a solution procedure for this class of problems by considering a common cycle length for all items produced on the machine. In this approach, in every production cycle, i.e., common cycle (T), each item was produced once. This method has a simple solution procedure and guarantees finding a feasible schedule. It was shown by Jones and Inmann (1989) that common cycle method performs well in many real world problems by providing optimal or near optimal solutions. Following these studies, Johnson and Montgomery (1974) proposed an exact solution procedure for multi-product single machine problems under assumption of common cycle when there is no priority for the production of items. Then, Haji et al. (2008) studied a production system consists of a machine producing defective items and considered non-zero setup times for rework cycle in the model. Afterwards, Taleizadeh et al. (2010) studied an EPQ problem for an imperfect production producing with a repair failure. Then, Taleizadeh et al. (2014) considered maintenance policies and modelled an imperfect production system with shortages for a single machine problem. In a more recent study, Pasandideh et al. (2015) investigated a multi-item economic production model for an imperfect production system with several types of failure. In this study the shortage was allowed in backorder form. Nobil and Taleizadeh (2016) proposed a multiitem single machine problem for an imperfect production system where the faulty items can either be repaired by rework or be sold in their current shape with a discounted price. Finally,  extended the single machine problem, proposing an EPQ problem for an imperfect production system consists of multi-machine with wastages, assuming some items can be produced on several machines.
Since chain of producing items, i.e., material flow from suppliers to final products delivery to customers, is one of the main issues in inventory control systems and supply chain problems, some researchers have developed models that consider the supply chain as a whole. For instance, Wee and Widyadana (2013) investigated a single-vendor singlebuyer inventory model with lost sales. Machine's unavailability time was considered to be two types; uniform and exponential distributions and orders were delivered in batches. They could derive some conditions which guaranteed the convexity of the model. Jha and Shanker (2013) studied a supply chain consists of a wholesaler that manufactures a product in batches and sell it to a number of customers. The demand rate of each customer followed an independent normal distribution and its lead time was reducible at a cost. The customers employed a continuous review policy, and the excess demand was fully backordered. The optimal combination of production-inventory policy of the system obtained by optimising collaborative expected cost of the system. Subsequently, Glock and Kim (2014) proposed a model in which a single buyer buys a product from multiple vendors with shipment consolidation. The study showed that consolidating shipments to the buyer from the vendor groups may lead to a significant reduction in total cost. Finally, they proposed a three-step solution procedure using steepest-descent approach in order to solve the proposed problem. Next, Muniappan et al. (2015) proposed optimal order quantity for a perishable production-inventory model in presence of rework and backordering under two different scenarios, coordination and incoordination. They assumed when coordination is in effect, the customer itself manages screening and disposal process, discount is available but shortage is not allowed. On the other hand, when there is no coordination in effect, producer is responsible for screening and disposal process but shortage is allowed for customers. Treviño-Garza et al. (2015) studied a vendor-buyer with defective items considering two different cases for lot size, discrete and continuous. The optimum solution was derived using discrete shipments value which is a more realistic assumption. In a recent study, Seifbarghy et al. (2016) proposed a two-level supply chain model consists of an item, a single producer and multiple retailers that were coordinated by a vendor-managed policy. Moreover, producer employed EPQ policy and some contracts based on the retailers' sales number, sales value and a set production rate. Due to the nature of the problem, a nonlinear integer programming model employing a discrete particle swarm algorithm was used to solve the problem. Next, Hariga et al. (2016) studied a wholesaler and a retailer coordination as a mixed-integer non-linear programming model, where packaging was performed employing returnable transport items (RTIs), for instance pallets, to guarantee its safety. After emptying the RTIs, they would be repaired and returned to the wholesaler. However, due to the unforeseen events, e.g., damages, the return time of the RTIs was assumed to be random and a rental option for providing them was considered.
Based on the literature, we can say that most studies in production planning and inventory control systems concentrated on vendor, i.e., manufacturer or wholesaler; and customer, i.e., retailer or seller; relationship and the raw material and consumers were neglected. In this study, we investigate an inventory control system which includes three levels, i.e., supplier, producer and distributor, with several raw materials and products. The raw materials shipped in a fixed size, e.g., pallets, trucks and etc., then produced on a single machine and sent to distributors in packages. We have developed two heuristics based on GA and biogeography-based optimisation (BBO) algorithm in MATLAB environment to both solve the problem and validate the accuracy of its generated result. Compared to other researches in the field, our approach is more practical rather than theoretical.

Problem definition and notation
As mentioned in the Literature review, Wee and Widyadana (2013) investigated a singlevendor single-buyer inventory model with lost sales, discrete delivery orders and stochastic machine unavailability time. Glock and Kim (2014) investigated a model in which a single buyer orders a product to multiple vendors. The vendors deliver their products to the buyer as joint shipments.
In this paper, we extend the two aforementioned models considering a multi-product three-echelon supply chain including multiple suppliers, single manufacturer and single distributor with discrete delivery orders for raw materials and products. The manufacturer produces n products through n types of raw materials, and sends them to the distributor. The required raw materials are replenished from the suppliers (i.e., product i requires raw material i which is received from supplier i). In this model, we assume that a part of the batches sent to the distributor are defective due to imperfect manufacturing systems, inefficient vehicles, improper packaging and etc. The defective products are scrapped after inspection by the distributor at the end of each shipment period of the distributor.
The notations and assumptions of the model are given as follows: The mathematical model is developed based on the following assumptions: 1 since all products are manufactured on a single machine with a limited capacity based on joint production strategy, an identical cycle time for all products is considered; i.e., 1 2 11 2 based on relationship between consumption of each raw material and the manufactured product, the cycle time of raw material i is equal to the corresponding production uptime which means 3 the first shipment of raw material i arrives when product i is to be produced (see Part A and B in Figure 1). Figure 1 illustrates the on-hand inventory of the system considering single-supplier (i.e., single raw material), the manufacturer and the distributor. Part A in Figure 1 represents the raw material inventory of the supplier where the supplier provides the raw material and sends to the manufacturer. Part B illustrates the raw material inventory of the manufacturer, wherein raw material is converted into a product with consumption rate of (-P); i.e., the system needs one unit of raw material to manufacture one unit product. Part C illustrates the product inventory of the manufacturer send to the distributor in batches. Part D represents the product inventory of the distributor, wherein the product is consumed with demand rate of (-D). Initially, based on Figure 1, the production uptime and downtime for product i is given in equations (1) and (2).

Average inventory analysis
(2) Figure 1 On-hand inventory of the system with single product The consumption cycle time of raw material i at the manufacturer is given as in equation (3).
The following equations can be obtained from Figure 1: where μ i represents the production quantity during the production time.
Based on Part A and B of Figure 1, the equation is obtained: where this equation is clearly shown in Figure 2. Replacing Tp i and Td i in equation (7) with the given values in equations (1) to (4), the equation (8) is obtained: The average inventory of raw material i at supplier i (ISR i ) can be given by equation (9).
( 1 0 ) The average raw material i inventory of the manufacturer in one cycle can be given as in equation (11).
( 1 2 ) Referring to Wang and Sarker (2006), the average product i inventory of the manufacturer in a single cycle is given as in equation (13).
( 1 4 ) The average product i inventory of the distributor in one cycle is given as in equation (15).

Total cost of the whole supply chain
The total cost function of the chain is the sum of the setup (CS), ordering (CO), producing (CP), transportation (CT), inspection (CI), scrapping (CD) and holding costs (CH). It can be computed from equation (17).

TC CS CO CP CT CI CD CH
The total setup cost of the suppliers and the manufacturer is given as in equation (18).
The total ordering cost of the suppliers and the manufacturer is given as in equation (19) The total production cost is the sum of production cost of all products and production cost of all raw materials. The total production cost is given as in equation (20).
The total transportation cost is the sum of the transportation costs of all products and all raw materials. It is given as in equation (22).
The total inspection cost is the sum of the inspection costs of all shipments. It is given as in equation (23). Using equation (5), equation (23) can be stated as in equation (24).
α i Q i is the quantity of scrap items in each delivered order. Thus, the total scrapping cost can be given as in equation (25).
In the joint production systems, the sum of production and setup times must be smaller than the cycle length. Which means, ∑ ∑ must be less than or equal to T. This constraint of the model is given as in equation (28).
Then, equation (28) can be stated as in equation (29). Through Constraint (31), we ensure that the delivery times of raw material i at supplier i is less than or equal to the production uptime of product i. Considering equation (8), constraint (32) can be stated as Constraint (33).
Based on the problem's assumptions, this constraint is redundant. In the joint production systems, the sum of consumption time of raw material and first shipment period of raw material i should be smaller than the cycle length. Which means, Tp i + Ts i Tp i must be less than or equal to T. This constraint is shown in (34).  (1) and (8), constraint (34) can be stated as Constraint (35).
Based on the definition of the problem, the number of shipments of product i during the cycle time must be larger than the number of shipments of product i during the production time; this can be given by (36). The number of shipments of raw material i during the time period must be smaller than or equal to delivery orders, given by (37). 3 7 ) The number of shipments of raw material i during the time period should be smaller than delivery orders, demonstrated given by (38). The last constraint specifies time period as in (40).The maximal value of T must be the time horizon.

Solution methods
In this section, the solution procedure for proposed MINLP model is presented. It contains large numbers of variables, constraints and a non-linear objective function. Such attributes make the problem complex and hard to solve using exact methods. Hence, we used GA to solve the problem, which is efficient and extensively deployed in the field of inventory problems (e.g., Pasandideh and Niaki, 2008;Pal et al., 2009;Moin et al., 2011;Taleizadeh et al., 2013;Saracoglu et al., 2014). We also used BBO to validate GA's results.
GA is a heuristics formed based on stochastic search and imitate the natural selection process. The most significant difference between GA and traditional search procedures is using a set of randomly generated solutions known as population. Each one of these solutions, which is called a chromosome, evolve to form a new generation. The chromosomes are assessed by the criteria of fitness function; the fitter chromosome has more chance to take part in breeding next generation. Moreover, to reproduce a new chromosome from its parents, there are two kinds of operators called crossover and mutation. Convergence of the algorithm happens after a number of iterations.
On the other hand, BBO is another evolutionary algorithm which searches the best solution using biogeography's principles: speciation, migration and extinction. MacArthur and Wilson (1967) proposed a mathematical model for this problem. Thereafter, Simon (2008) introduced an adaptation of BBO logic for optimisation problems. In this method, the solutions are represented by islands; these islands have a population and immigration rate. The solutions improve by either immigrating or mutation operator.

Proposed GA-based heuristic
The proposed GA-based heuristic's steps for solving the MINLP are as follows.

Decision variables' boundaries
The first step in solving the problem is to obtain the upper and lower bounds of decision variables. These boundaries reduce the solution space and subsequently lessen the computational time. These boundaries can be obtained as follows: • Boundaries of T: based on (49), the acceptable interval for this variable is [0, 1].
• Boundaries of x i : considering constraint (48) the lower bound of x i is equal to 1. On the other hand, based on constraint (46), we have: According to constraint (50) the upper bound of x i is represented by S i J .
• Boundaries of w i : with respect to constraints (45) and (48), w i is higher than or equal to 2. On the other hand based on constraint (47), we have: According to constraint (51) the upper bound of w i is represented by P i J .
• Boundaries of m i : based on constraints (45) and (48), m i is higher than or equal to 1; On the other hand, based on constraints (45) and (51) we have:

Initial definition
Chromosome: every solution to the proposed MINLP is called a chromosome.
• Generation population: equals to the number of chromosomes showed by nPoP.
• Crossover probability (Pc): shows the probability of crossover occurring for a chromosome. Generally, in each generation the number of crossovers can be obtained by Pc × nPoP.
• Mutation probability (Pm): shows the probability of mutation occurring for a chromosome. Generally, in each generation the number of mutations can be obtained by Pm × nPoP.

Chromosome representation
The proposed chromosome contains four decision variables including T, x i , m i , w i . As a result, this chromosome has 3n + 1 genes, in which one gene used for T and n genes for each one of the remaining variables. An example of this chromosome illustrated in Figure 3. The initial population, including nPoP individuals, is randomly generated in two steps. The first step consists of generating a random number (=k) from [0, 1] interval for genes of these chromosomes. In the second step, these numbers are normalised by adding the lower bounds of each gene to a product of gene's interval length and k.

Penalty function
In this algorithm the chromosomes with less fitness function values are more favourable.
Since the constraints in inequalities (42) to (45) could adversely affect the solution's feasibility, we add a penalty to objective function in terms of potential violations. This penalty performs based on equation (53), where M, g(x), and P(x) represent a reasonably large number, the constraint under investigation (g(x) ≤ b), and penalty function of infeasible solution (chromosome), respectively. Afterwards, updating cost function is done using equation (54), where A(x) represents the fitness function value of a solution.

Crossover
Note that, in this algorithm the single point crossover operator is used to generate new chromosomes. For the first step, two chromosomes (parents) are randomly selected from former generation then employing crossover operator two new chromosomes (offspring) are obtained. To obtain offspring, T part of parents are swapped then three random numbers from {1, 2, …, n -1} will be chosen to obtain cutting points of x, m, w parts of chromosomes. Finally, we will exchange these parts of chromosomes to form offspring; this procedure is performed in a generation to achieve the Pc × nPoP new chromosomes. An example of the mentioned procedure is illustrated in Figure 4.

Mutation
Initially Pm × nPoP chromosomes are chosen randomly from last generation then using random mutation, new chromosomes will be obtained as follows: • Ω is a randomly generated number from 1 to 4 (the number of decision variables).
• Choose Ω decision variables by chance.
• Based on the selected genes, some genes are chosen and reproduced randomly. An example of the mentioned mutation operator demonstrated in Figure 5.

Stopping criteria
After minIt iterations if the algorithm converged, it will stop. The algorithm must be run minIt times to prevent early convergence.

Proposed BBO-based heuristic
BBO inspired from natural investigations of species distribution during different eras among diverse places (Simon, 2008(Simon, , 2009). This population-based algorithm uses a set of candidate solutions to solve an optimisation problem. In this algorithm each habitat represents a candidate solution (an individual) and has its habitat suitability index (HSI) to show the degree of suitability. In other words, higher-HSI habitat characterises better solutions and vice-versa. Solution features emigrate from high-HSI habitats (emigrating habitat) to low-HSI habitats (immigrating habitat), however this emigration is not equivalent to their distinction. This is because low-HSI habitats have the potential to accept a large number of new features from high-HSI habitats through an immigration process. Consequently, emigration and immigration operators are employed to improve solutions of this optimisation problem.
Just to illustrate, let us consider Figure 6 in which curve I and E denote the maximum possible immigration and emigration rates, respectively. In this setting, S max indicates the largest possible numbers of species that a habitat can support (where the immigration rate equals to zero) while S 0 is the point where immigration and emigration rates are equal.
In this study a BBO is employed to assess performance of the proposed GA. In spite of the tendency to consider similar operators for both BBO and GA, there are some differences due to their features. A brief comparison between proposed GA and BBO is as follows. Initialisation phase of proposed BBO consists of random distribution of species among diverse habitats (=chromosomes) with specific characteristics (=genes) and calculating overall rating of suitability index (=fitness function). In the BBO the population of a habitat balances based on a migration operator. This operator acts as a function of immigration and emigration probabilities, presented by λ and μ respectively.
These probabilities are a function of the habitat population, for a habitat with a great deal of species due to the surviving chance reduction and resource limitations, immigration and emigration probabilities will decrease and increase, respectively (see S 1 and S 2 in Figure 6). On the other hand cataclysmic operator can change the habitats' suitability indices. Thus it can be concluded that the migration and cataclysmic operators act similar to the crossover and mutation operators.

Computational results
This section provides an analytical framework for the proposed GA and BBO algorithms and represents numerical results. In the subsection, the input parameters of the MINLP are chosen randomly from Table 1. Table 1 Initial parameters of the MINLP problem

Parameters tuning
In this study, a number of normalised responses such as fitness function (Y 1 ) and CPU time (Y 2 ) are optimised using response surface methodology (RSM) with regard to initialisation parameters, for instance, population size (X 1 ), crossover probability (X 2 ), mutation probability (X 3 ), and minimum number of iteration (X 4 ).  In the beginning, experiments were performed on 20 treatment combinations based on the levelled variables in Table 2, where every individual response was the average result of ten times running MINLP with six items (see Tables 3 and 4). In addition, the time horizon is considered one year for this problem. Afterwards, fitting the responses were performed using corresponding parameters according to equations (55) and (56). Moreover, employing ANOVA the accuracy of fitness functions were examined (Tables 5 and 6). At the end, related parameters for suggested GA were derived by solving optimisation model of equation (57) (see Table 7). Besides, BBO calibrated parameters can be determined using the same procedure (shown in Table 8).  5 7 ) where W j denotes the weights of the j th objective in which W 1 + W 2 =1. We assumed W 1 = 0.8and W 2 = 0.2.

Numerical results
The proposed GA and BBO algorithms were run on a personal computer equipped with 2.4 GHz CPU and, 4 GB of memory and MATLAB 2014a software. These algorithms were run 10 times on 10 different problems, the average of CPU time and total costs are shown in Table 9. The input parameters of these algorithms were similar to the ones shown on Table 1. The input parameters were randomly selected from Table 1 and incorporated into the algorithm parameters shown in Tables 7 and 8. Based on the total cost and time, it can be concluded that GA performed slightly better than BBO. Although non-parametric Kruskal-Wallace test showed no significant difference between the two algorithms (see Table 10).

Table 10
Summary of Kruskal-Wallis test for two measures according to information given in Table 9 Measure P-value Result

Sensitivity analysis
This subsection shows result of the sensitivity analysis performed on important parameters of the problem (i.e., raw materials' production rate, manufacturing rate of end product, demand's rate and rate of failure). In the process, each parameter was changed while maintaining other parameters' original value (see Table 11). As it can be seen from Table 11, the proposed MINLP has the least sensitivity to R i .

Conclusions and future research
In this study, a non-linear three-echelon inventory model was considered in which the delivery of ordered raw materials and final products were done in batches. The addressed system in this paper was a three echelon supply chain including a supplier, a producer and a distributor. Moreover, based on real world instances that items are shipped in discrete sizes, we assumed the raw materials and products are shipped in batches and the objective was to find optimal batch size and shipping time so as to minimise the total costs. Each final product needs a specific raw material which is produced by a specific supplier. The studied model was under single machine assumption with a limited joint production capacity for the cycle length of all products in the manufacturer. Considering our practical assumption of imperfect production system, the final product needs inspection in order to separate defective items and enhance customer satisfaction. The number of shipments of each raw material from the supplier to the manufacturer, the numbers of shipments of each product from the manufacturer to the distributor and the cycle length were determined in such a way that minimises the inventory system, including setup, ordering, producing, transportation, inspection, scrapping and holding costs. We adopted a more practical approach by developing two heuristics based on GA and BBO algorithm in MATLAB environment to both solve the problem and validate the accuracy of its generated result. The GA-based heuristic was used to solve this problem and the BBO-based heuristic was employed to evaluate its validity. Moreover, RSM was utilised for tuning the parameters of the heuristics. Finally, a slightly better solution for the proposed GA was obtained compared with the proposed BBO.
For further studies, the following developments are proposed: • considering multiple suppliers for each raw material • considering multiple distributors in the model

• considering storage constraint
• studying the model when reproduction is allowed • studying the model when shortage is allowed • considering maintenance policies for machine • studying a relatively same system when rework is permitted • considering marketing and pricing to enhance customer demands • model the system when postponement policy is allowed by producing large amounts of items in beginning of the period and postponing some to the end of the period • each machine has its own unavailability times.