- Email: [email protected]
A genetic algorithm for a joint replenishment problem with resource and shipment constraints and defective items
Int. J. Production Economics 175 (2016) 142–152
Contents lists available at ScienceDirect
Int. J. Production Economics journal homepage: www.elsevier.com/locate/ijpe
A genetic algorithm for a joint replenishment problem with resource and shipment constraints and defective items P. Ongkunaruk a, M.I.M. Wahab b,n, Y. Chen c a
Department of Agro-Industrial Technology, Kasetsart University, 50 Ngam wong wan Road, Lad Yao, Chatuchak, Bangkok 10900, Thailand Department of Mechanical and Industrial Engineering, Ryerson University, 350 Victoria Street, Toronto, ON, Canada 5MB 2K3 c School of Economics and Management, Southwest Jiaotong University, Chengdu, Sichuan 610031, PR China b
art ic l e i nf o
a b s t r a c t
Article history: Received 3 June 2015 Accepted 18 February 2016 Available online 2 March 2016
A joint replenishment problem (JRP) is presented to determine the optimal reordering policy for multiple items with defective quantity and several restrictions such as shipment constraint, budget constraint, and transportation capacity constraint. The objective is to minimize the total expected cost per unit time. A two-dimensional genetic algorithm (GA) is provided to determine an optimal family cycle length and the reorder frequencies. A numerical example is presented and the results are discussed including the effect of defective items on the ordering policy. Extensive computational experiments are performed to test the performance of the GA. The JRP was also solved by a differential evolution (DE) algorithm and the results obtained from DE were compared with those obtained from GA. Crown Copyright & 2016 Published by Elsevier B.V. All rights reserved.
Keywords: Joint replenishment Defective items Shipment constraint Genetic algorithm Differential evolution algorithm
1. Introduction A joint replenishment problem (JRP) is defined as the coordination of the replenishment of a group of items that may be ordered jointly from a single supplier (Goyal, 1974). One of the advantages of jointly placing order for a group of items is sharing the fixed cost associated with order, and hence the total cost is reduced when compared to the cost of individually placing order for each item (Moon and Cha, 2006). Moreover, a JRP is a good practice to reduce transportation cost, when a vehicle is used to transport multiple items simultaneously. Similarly, the supplier's set-up cost can be reduced by manufacturing multiple items using the same facility. One of the disadvantages of the JRP is increasing the average inventory level and system control costs. In order to solve the basic JRP, family cycle time and the frequency of replenishing individual items have to be determined while minimizing the total cost, which is the sum of ordering and holding costs. Most of the studies related to the basic JRP have usual assumptions such as deterministic and constant demand, linear holding cost, no shortage, and no quantity discount (e.g., Khouja and Goyal, 2008; Robinson et al., 2009). The JRP has received considerable attention over the last few decades and different solution approaches have been proposed to solve the JRP. For example, Shu (1971), Goyal (1974), and Silver n
Corresponding author: Tel.: þ 1 416 979 5000x2670; fax: þ1 416 979 5265. E-mail addresses: [email protected] (P. Ongkunaruk), [email protected] (M.I.M. Wahab), [email protected] (Y. Chen). http://dx.doi.org/10.1016/j.ijpe.2016.02.012 0925-5273/Crown Copyright & 2016 Published by Elsevier B.V. All rights reserved.
(1975) proposed a simple procedure to determine the order quantity for multiple items in a JRP. Moreover, many heuristics have been developed to solve different versions of JRP, namely, basic JRP, JRP with quantity discount, and constrained JRP using various methods including direct and indirect grouping strategies, genetic algorithm (GA), evolutionary algorithm (EA), RAND, QD– RAND, and C-RAND (e.g., Kaspi and Rosenblatt, 1983, 1991; van Eijis et al., 1992; Hariga, 1994; Wildeman et al., 1997; Khouja et al., 2000; Cha and Moon, 2005; Moon and Cha, 2006; Olsen, 2008; Hong and Kim, 2009). However, among many heuristics, GA and EA have been proved to be effective for solving the JRP (Olsen, 2008). Moreover, it was concluded that solutions obtained from GA combined with a local search algorithm and RAND are very close to the global optimum (Hong and Kim, 2009). Recently, Wang et al. (2012) proposed a differential evolution (DE) algorithm based on Storn and Prince (1997) to solve a JRP and claimed that the DE algorithm can solve a JRP more effectively than the EA used in Olsen (2008). Hence, in this research, we used GA and DE, which are the best algorithms existing in the literature, to solve a JRP. With respect to GA literature, we used a two dimensional chromosome representation. There are several studies on twodimensional GA for solving packing problem (e.g., Herbert and Dowsland, 1996; Bortfelldt, 2006), cutting problem (e.g., Ono and Ikeda, 2008), and knapsack problem (e.g., Bortfelldt and Winter, 2009). In addition, the two-dimensional GA was also implemented in different areas of research, for example, ising problem (Anderson et al., 1991) and photonic crystal slab (Jia et al., 2009). With respect to DE, Storn and Prince (1997) proposed a novel DE to
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
solve a complex continuous non-linear functions and we used a similar algorithm to solve the JRP. Most of the studies on JRP assumed that items included in an order can be shipped together without any restrictions. However, in practice, there are restrictions or dependency among the items. As a result, though an order is placed for certain items at the same time, they cannot be shipped together in the same truck. For example, chemical or fragrant products should not be shipped together with dry food or food ingredients, since the odor will be absorbed by the dry food or food ingredients. Another example is ambient products, which are not shipped with fresh produce and ready to eat meals due to different storage temperatures. Recently, several studies (e.g., Olsen, 2008; Wang et al., 2012) considered such interdependency of the items in a JRP. However, these studies assumed that such interdependency increases minor ordering costs to maintain different requirements of the prohibited items in the same truck and then the JRP was solved by using direct and indirect groupings. Let us consider a modern retail business company, which operates in Thailand, having restrictions in its shipment policy to maintain quality and food safety. This company does not allow transporting prohibited items in the same truck at all. This retail company plans its ordering policy to minimize its total inventory related cost considering quality and safety requirements of the company. Then, the company provides the ordering policy to a third-party logistics provider and it ships the items for a fixed transportation cost. Hence, the JRP models in Olsen (2008) and Wang et al. (2012) cannot be applied to this retailer. In order to solve the retailer's JRP, in this paper, without changing minor setup costs of items having restrictions, shipment restrictions or interdependency of the items are included in the constraints and then a constrained JRP is solved. In addition, there are several other constraints besides shipping restrictions or interdependency of items. For example, Goyal (1975) introduced a JRP with budget constraint. Khouja et al. (2000) addressed the issue of maximum load constraint in a JRP. Later, Moon and Cha (2006) also studied a JRP with a budget constraint. However, the existing literature on JRP does not consider such constraints simultaneously in a JRP. In this research, we included budget constraint and transportation capacity constraint with interdependency of the items and then solved a constrained JRP. Another important aspect in a JRP is defective items. Wahab et al. (2011) pointed out that many Thai small and medium enterprises implemented total quality management, however buyers still sometimes find defective items in the shipments, and it is also possible that some items may also get damaged during shipment. There are several studies dealing with defective items in a single item replenishment such as EOQ, for example, Porteus (1986), Rosenblatt and Lee (1986), Salameh and Jaber (2000), Cárdenas-Barrón (2000), Goyal and Cárdenas-Barrón (2002), Papachristos and Konstantaras (2006), Maddah and Jaber (2008), and Wahab and Jaber (2010). The defective items in shipments will affect the ordering policy. For example, if there are more defective items in the shipment, the order quantity should be higher to compensate for those defective items. Hence, this study also considered the percentage of defective items in the JRP to investigate its effect in the ordering policy. The purpose of this research is to model a constrained JRP for multiple items, where each item has a certain percentage of defective item in each shipment, with constraints such as dependency among items (i.e., some items have restrictions that they cannot be shipped together in the same truck), budget constraint, and transportation capacity constraint. Such a JRP exists in many real business applications. For example, this problem can be found in a retailer that delivers products from its distribution center to retail stores. The JRP is solved by using a two-dimensional GA and DE. The results obtained from both GA and DE are compared.
143
This paper is organized as follows: Section 2 presents the model, notations, and decision variables. Section 3 presents the solution methodology, a two-dimensional GA, and its steps. Section 4 illustrates numerical examples and also presents the comparison between GA and DE. Section 5 concludes the paper.
2. The model This section first presents indices, parameters, and decision variables of the model and then describes the JRP step by step to develop the model. 2.1. Indices, parameters, and decision variables Indices i an index of an item, where i ¼ 1; 2; …; n j an index of a truck that is used to ship items, where j ¼ 1; 2; …; J l an index of a set representing items that can not be shipped together, where l ¼ 1; 2; …; L Parameters n the number of items J the number of trucks L the number of sets of items that cannot be shipped together Δj the capacity of truck j Di Qi ai
the the the the
demand rate for item i order quantity for item i handling cost for item i unit screening cost for item i
μi νi the lost per defective unit for item i pi the percentage of defective of item i xi the screening rate for item i ti the screening time for item i hi the holding cost per unit of item i per unit time ci the unit variable cost for item i K the ordering cost per order B the amount of capital that can be invested Sl the lth set of items that cannot be shipped together Decision variables T time interval between replenishment for a set of items mi the integer number of intervals that the replenishment takes place for item i bij the proportion of item i transported by truck j ( 1 if item i is ordered and delivered by truck j; Yij 0; otherwise:
2.2. The joint replenishment problem Item i has demand Di per year. It is assumed that demand for each item is deterministic and constant, and shortages are not allowed. Item i is ordered in every mith replenishment cycles, where mi is a positive integer, meaning the replenishment quantity of item i will last for miT period of time. This also indicates that the cycle time for item i is miT. An important aspect that has been addressed in this problem is the restriction in shipping. That is, because of the nature of the products and requirements, some items cannot be shipped together. For example, chemical or fragrant products should not be shipped in the same truck as dry food or raw material. Though
144
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
these products were ordered at the same time, they are shipped in different trucks with other products that do not have any conflict. In other words, several items can be shipped simultaneously by a truck if there is no restrictions in shipping them together. In this scenario, depending on capacity requirements for shipping, sometimes an ordered quantity of an item may not be shipped by a single truck and it can be divided and shipped by more than one truck at the same time. A joint ordering cost K is incurred every T period of time. For example, the cost of placing an order for multiple items via fax, email, or phone, and this cost is independent of types of items and their order quantities. The handling cost ai incurs after receiving item i for every miT period of time. Therefore, the total setup cost per unit time is Kþ
Pn
i¼1
PJ j¼1
ai bij Y ij mi
T
;
n X
c i Di : 1 E½p i i¼1
Q i ð1 E½pi Þ Z Di t i
ð9Þ
The total lost of defective units, due to markdown, is given by n X νi Di E½p i
i¼1
1 E½pi
:
ð10Þ
The total expected cost per unit time for the JRP is the summation of all aforementioned costs and it is given as follows: 0
E½CðT; m0i s; bij s; Y 0ij sÞ ¼
n X Di ðci þ μi þ νi E½pi Þ 1 E½pi i¼1
Di : xi
ð2Þ
ð3Þ
ð5Þ
ð6Þ
By substituting for Qi and ti, the above equation can be expressed as: Di mi T D2 m TE½pi þ i i : 2 xi ð1 E½pi Þ2
þ
PJ j¼1
ai bij Y ij mi
T " # D i mi T D2 m TE½pi þ i i þ hi : 2 xi ð1 E½pi Þ2 i¼1
ð11Þ
This cost function needs to be minimized subject to capacity, budget, and shipping constraints as follows: n X Di mi bij Y ij T i¼1
ð1 E½pi Þ
r Δj ;
8j
ð12Þ
n X D i mi c i T r B; ð1 E½pi Þ i¼1
ð13Þ
bij r Y ij ;
ð14Þ
8 i; j
0 r bij r 1;
8 i; j
ð15Þ
bij ¼ 1;
8i
ð16Þ
Y ij r 1;
8 j; l
ð17Þ
J X
X i A Sl
T Z0
ð18Þ
mi A f1; 2; …g; Y ij A f0; 1g;
ð4Þ
where Qi consists of good and defective items, and their average inventory level is given by Q i ð1 E½pi Þ Q i E½pi t i þ 2 mi T
i¼1
j¼1
Since only good items can be used to meet the demand during miT periods of time, the replenishment quantity Qi for item i is given by Q i ð1 E½pi Þ ¼ Di mi T;
Pn
n X
Substituting t i ¼ Q i =xi , it becomes
Ii ¼
μi Di : 1 E½pi i¼1
ð1Þ
Moreover, the following condition must be satisfied for no-shortage:
Ii ¼
n X
Kþ
in this expression, the second term indicates that the handling cost is charged for each fraction bij of the order quantity shipped in each vehicle. In most industrial applications, due to the damage in transit and/or imperfect production process, some defective items present in a shipment. It is assumed that pi is the percentage of defective for item i in each replenishment, where pi is a random variable. Once the replenishment is received, it is screened at the rate of xi, which is greater than Di, during the time t i ¼ Q i =xi . When the defective items are segregated by the screening process, those defective items are immediately sold for a markdown price, and hence defective units of item i are retired from the system at time ti. For item i, the expected number of good units in each replenishment is Q i ð1 E½pi Þ, and only good units can be used to meet the demand in each cycle. Hence, the total purchasing cost per unit time is
E½pi r1
Since μi is the screening cost per unit, the total screening cost per unit time is given by
ð7Þ
A holding cost of hi for item i is incurred per unit time. Therefore, the total expected holding cost per unit time is " # n X DmT D2 m TE½pi hi i i þ i i : ð8Þ 2 xi ð1 E½pi Þ2 i¼1
8i
ð19Þ
8 i; j
ð20Þ
Eq. (12) is the capacity constraint. The total shipment size in a truck cannot be more than its capacity. Eq. (13) is the budget constraint. The total purchasing per order should not be more than the budget. Inequality (14) is to ensure that if Yij is zero, then bij is zero. Inequalities in (15) are the boundary constraints of bij. Eq. (16) implies that item i can be fully or partially delivered by truck j. Inequality (17) is the shipment requirement constraint. Items in set Sl cannot be shipped simultaneously in the same truck. Eq. (18) implies that the time interval must be positive. Eq. (19) implies that mi must be positive integer that is greater than or equal to one. Eq. (20) implies that Yij is binary. The objective is to determine the decision variables, mi, T, bij, and Yij by minimizing the total expected cost per unit time given in Eq. (11). Therefore, taking the first derivative of the total expected cost with respect to family cycle length T and setting it to zero, we obtain Kþ
0
dE½CðT; m0i s; bij s; Y 0ij sÞ dT
Pn
¼
i¼1
PJ j¼1
ai bij Y ij mi
T2 # Di mi D2i mi E½pi þ þ hi ¼ 0: 2 xi ð1 E½pi Þ2 i¼1 n X
"
ð21Þ
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
145
After rearranging the above equation, we obtain the following solution for T and it is denoted as T0. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P P ai bij Y ij u K þ ni¼ 1 Jj ¼ 1 u u mi " #: T0 ¼ u ð22Þ 2 uP t n h Di mi þ Di mi E½pi i¼1 i 2 xi ð1 E½pi Þ2 From the capacity constraint, Eq. (12), it can be written as: T 1 ¼ min
j ¼ 1;…;J Pn i¼1
Δj Di mi bij Y ij ð1 E½pi Þ
:
ð23Þ
From the budget constraint, Eq. (13), it can also be written as: T2 ¼
B : Di mi ci i¼1 ð1 E½pi Þ
ð24Þ
Pn
0
Finally, for a given set of m0i s, bij s, and Y 0ij s, the optimal cycle time Tn can be determined as follows: T n ¼ minfT 0 ; T 1 ; T 2 g:
ð25Þ 0 bij s,
Y 0ij s
and that The next step is to determine the best set of m0i s, minimize the total expected cost per unit time given in Eq. (11). In 0 order to determine m0i s, bij s, and Y 0ij s, GA is employed in the next section.
3. Genetic algorithm Genetic algorithm starts with an initial set of random solutions called a population. Each solution in the population is called a chromosome. Each chromosome composes a set of genes. All genes are generated by the random number generator of uniformly distributed random number between 0 and 1. For each generation, the chromosomes are evaluated using some measures of fitness. The algorithm is presented in Fig. 1 and explained as follows: 3.1. GA chromosome representation Fig. 1. GA.
The basic cycle time T, n integer values of m0i s, n J values of and n J values of Y 0ij s have to be decided for solving the JRP. In order to determine the values of those decision variables, the following steps are given with examples. 0 bij s,
1. Encode the chromosome that has n ðJ þ 1Þ genes. Each gene represents a uniform random variable between 0 and 1. An example is given in Fig. 2 for i¼ 5 and j¼ 2. The first column represents the values of m0i s and the rest represent the values of 0 bij s. 2. Decode the chromosome that represents the solution of the JRP as follows: GeneðiÞð1Þ that represents mi can be decoded in the following way: UB LB mi ¼ mLB i þ⌊ðmi mi þ1Þ GeneðiÞð1Þc UB
8 i ¼ 1; 2; …; n;
ð26Þ
LB
where the upper bound, mi , and a the lower bound, mi , of mi can be determined from Eqs. (33) and (34), respectively. Then, there are two cases to be considered as follows: (a) If there is no restriction that items can be shipped in the same 0 truck, GeneðiÞð2Þ with n ðJ þ 1Þ genes that represent bij s can be decoded in the following way: bi1 ¼ GeneðiÞð2Þ
8 i ¼ 1; 2; …; n:
ð27Þ
Fig. 2. Encoding chromosome with n ðJ þ 1Þ genes.
bij ¼
1
j1 X
! bik
k¼1
GeneðiÞðj þ1Þ
biJ ¼ 1
J1 X
bij
8 i ¼ 1; 2; …; n; j ¼ 2; …; J 1:
8 i ¼ 1; 2; …; n:
ð28Þ
ð29Þ
j¼1
Next, determine Y 0ij s if bij 4 0; then Y ij ¼ 1; otherwise, Y ij ¼ 0. An example is given in Fig. 3.
146
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
(b) If there is restriction that some items cannot be shipped together in the same truck, then Y 0ij s, which do not violate Eq. (17), are randomly determined. Next, Gene ðiÞð2Þ with n 0 ðJ þ1Þ genes that represent bij s can be decoded in the following way: bi1 ¼ GeneðiÞð2Þ Y i1
bij ¼
1
j1 X
8 i ¼ 1; 2; …; n:
GeneðiÞðjþ 1Þ
k¼1
Y ij
biJ ¼ 1
8 i ¼ 1; 2; …; n; j ¼ 2; …; J 1:
J1 X
bij
8 i ¼ 1; 2; …; n:
ð31Þ
ð32Þ
j¼1
An example is given in Fig. 4. It is well-known that if the cost is minimum at mi, then it is lower than the cost at mi 1 or at mi þ 1 . Hence, in order to minimize the cost function given in Eq. (11), the lower bound of mi, miLB, and the upper bound of mi, miUB, can be defined from the optimality condition as follows (e.g., Goyal, 1974): "
UB mUB i ðmi 1Þ r
hi T 2min
ai D2i E½pi
Di þ 2 xi ð1 E½pi Þ2
UB # r mUB i ðmi þ 1Þ;
8 i;
ð33Þ "
LB mLB i ðmi 1Þ r
hi T 2max
ai Di D2i E½pi þ 2 xi ð1 E½pi Þ2
i
ð36Þ
i
ð30Þ
! bik
9 8 > > > > > > >vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi> > > = > i ¼ 1;…;n>u > Di E½pi > > t Di > > > ; : hi 2 þ x ð1 E½p Þ2 >
LB # r mLB i ðmi þ 1Þ;
8 i;
ð34Þ where Tmax can be estimated from Eq. (22) by substituting mi equal to 1; Tmin can be obtained from the standard EOQ for an individual item. They can be expressed as given below: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P P u K þ ni¼ 1 Jj ¼ 1 ai bij " #: ð35Þ T max ¼ u u uPn Di D2i E½pi t þ i ¼ 1 hi 2 xi ð1 E½pi Þ2
3. Determine the optimal basic cycle time, Tn, from Eq. (25) by substituting a given set of m1 ; m2 ; …; mn , b11 ; …; bnJ , and Y 11 ; …; Y nJ . 4. Check the feasibility of the solutions by substituting them in Eqs. (12), (13), and (17). 5. Evaluate the fitness function. Calculate the total expected cost per unit time by using the objective function given in Eq. (11). Then, the fitness function is the summation of the objective function and penalty cost. There are three constraints: vehicle capacity, budget, and dependency of items. If the chromosome violates any constraint, then a larger penalty cost is added to the fitness function. Hence, the probability of being selected to be a parent is reduced. Otherwise, no penalty cost is added to the fitness function. 3.2. Selection We implement GA to select the best method by applying two selection methods: (i) a combination of rank and roulette wheel selection (Kruse et al., 2013), and (ii) an unbiased tournament. (i) Using the first method, initially the fitness function is ranked in the ascending order, and then the probability of selection is calculated. The probability of selection is determined by dividing its rank by the total summation of the rank. (ii) Unbiased tournament randomly divides the chromosomes into groups and selects the best solution from each group. 3.3. GA operators GA operators create a new generation from parent chromosomes. There are three GA generators as follows: (i) Elitism: In order to maintain the best solutions from the last generation to the next generation, a certain number of best chromosomes will be copied to the next generation with a specific probability (Beasley et al., 1993; Kruse et al., 2013). (ii) Crossover: The crossover of the selected parents will produce the child for the next generation. We apply horizontal crossover by randomly selecting any row of two parent chromosomes and exchange the rows (e.g., Anderson et al., 1991). An example is given in Fig. 5.
Fig. 3. Chromosome decoding when there is no restriction on the delivery of items.
Fig. 4. Chromosome decoding when item 1 can not be shipped with item 2, and item 4 can not be shipped with item 5.
Fig. 5. An example of horizontal crossover.
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
147
Table 1 Values of the parameters.
Fig. 6. An example of horizontal mutation.
(iii) Mutation: In order to prevent the local optimum, the offspring from crossover will be mutated with a certain probability. We apply horizontal mutation where a row of the chromosome will randomly be selected, and then two genes will be swapped randomly (Beasley et al., 1993; Kruse et al., 2013). An example is given in Fig. 6.
Qi ¼
D i mi T ; ð1 E½pi Þ
8 i ¼ 1; …; n:
ð37Þ
Note that when substituting pi ¼ 0; 8 i, the case with no defective items can be analyzed.
4. Numerical examples In order to illustrate the developed model numerically, a JRP is considered with 5 items in which item 1 must not be shipped with item 2, and item 4 must not be shipped with item 5; and there are two trucks each with capacity of 5000 items. The available budget, B, is $40,000 and the ordering cost, K, is $10 per order. The percentage of defective of item i, pi, is assumed to be uniformly distributed, i.e., pi U½vi ; wi , hence E½pi ¼ ðvi þ wi Þ=2 and Var½pi ¼ ðwi vi Þ2 =12. Moreover, the values of the other parameters are given in Table 1: The problem is solved by GA and the results are summarized in Table 2. Items 2 and 5 are assigned to truck 1, items 1 and 4 are assigned to truck 2, 35% of the order quantity of item 3 is assigned to truck 1 and the rest is assigned to truck 2. The constraint that item 1 must not be shipped with item 2, and item 4 must not be shipped with item 5 has been satisfied. Items 1 and 3 are ordered
1
2
3
4
5
Di (units/year) ci ($/unit) μi ($/unit) νi ($/unit) ai ($/delivery) hi ($/unit/year) xi (units/year) E½pi
1000 6 0.15 3 30 1.2 1500 0.01
4000 2 0.1 1 20 0.8 6500 0.02
3000 5 0.15 2.5 80 1 5000 0.015
5000 8 0.2 4 40 3.2 7000 0.04
2500 4 0.15 2 20 0.8 4000 0.030
Table 2 Summary of the results.
3.4. Initialization The initialization process randomly generates the first generation of chromosome equal to the number of population. The generated chromosomes will be decoded so that the solutions and fitness functions are determined. Next, the stopping criterion is checked. If it is not satisfied, then the new generation will be generated by two branches. During this process, first, the GA will preserve some good solutions through elitism by cloning the best chromosomes from the previous generation with a certain probability ðERÞ. Second, the parent chromosomes will be selected by two methods: RRW and TS. Third, the crossover will be performed to the selected parent chromosomes and the mutation will be applied ðMRÞ. This implies that some children after crossover may not be mutated, and finally, the total number of offspring will be equal to the number of population ðPÞ. Hence, at the combined arrow into the “new generation” in Fig. 1, P½MRð1 ERÞ þ ð1 MRÞð 1 ERÞ þ ER ¼ P. After that, the generation number is increased by one. The process will be repeated until the stopping criterion is satisfied as shown in Fig. 1. The stopping criterion is that the difference of fitness function between the current and the last generation is less than 5. The optimal value of family cycle time, T, the reorder frequencies, m0i s, the proportion of each item transported by each truck, bij, are obtained. The order quantity of item i is given by the following equation:
Item
Item
1
2
3
4
5
mi bi1 bi2 Yi1 Yi2
3 0 1 0 1
1 1 0 1 0
3 0.35 0.65 1 1
1 0 1 0 1
2 1 0 1 0
every 3T n , items 2 and 4 are ordered every T n , and item 5 is ordered every 2T n , where T n ¼ 0:079146 year. The expected total cost per year is $88,035.30. Now, in order to understand how defective items affect the ordering policy, we change the percentage of defective items for each item. It is assumed that vi ¼ 0; 8 i ¼ 1; 2; …; 5, and wi takes different values for different items, where 0:02 r w1 r 0:26, 0:05 r w2 r 0:29, 0:035 r w3 r 0:275, 0:03 r w4 r0:27, and 0:022 r w5 r0:262. It should be noticed that as wi is changed in Table 3, both mean and standard deviation of the uniform distribution are changed. From the first row of Table 3, it can be seen that when the items are perfect (i.e., vi ¼ wi ¼ 0), the order quantity for each item and the total expected cost are the lowest. As the percentage of defective items increases, the optimal family cycle length, T, becomes shorter indicating more frequent order. Moreover, order quantity of each item increases while mi for each item remains the same; and also when percentage of defective items increases the expected total cost per unit time increases. This implies that, while satisfying the constant demand for each item, in order to compensate for the defective items, the order quantity should be increased and also the items should be ordered more frequently. This leads to a higher cost when compared to the JRP with perfect items. Next, in order to evaluate the performance of GA, besides the total cost, three criteria are used:
Convergence score ðCSÞ: It determines how well the solution converges to the best value for all generations of GA (De Jong, 1980). CS can be computed by using the following equation: CS ¼
1 Gen
PGen
g¼1
E½C g ð .,.,.,. Þ
ming E½C g ð .,.,.,. Þ
;
ð38Þ
where Gen is the number of generations, E½C g ð .,.,.,. Þ obtained from Eq. (11) is the best total expected cost of generation g. CS C 1 means that the average fitness value is not different from the best value that occurred in this generation. Then it implies that there is no sufficient gradient to reach to the global optimum. Improvement score ðISÞ: This measure determines how well the solution was improved since the first generation (Beasley
148
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
Table 3 Sensitivity analysis for items with imperfect quality. w1
w2
w3
w4
w5
m1
m2
m3
m4
m5
T
Q1
Q2
Q3
Q4
Q5
E½C
0 0.020 0.040 0.060 0.080 0.100 0.120 0.140 0.160 0.180 0.200 0.220 0.240 0.260
0 0.050 0.070 0.090 0.110 0.130 0.150 0.170 0.190 0.210 0.230 0.250 0.270 0.290
0 0.035 0.055 0.075 0.095 0.115 0.135 0.155 0.175 0.195 0.215 0.235 0.255 0.275
0 0.030 0.050 0.070 0.090 0.110 0.130 0.150 0.170 0.190 0.210 0.230 0.250 0.270
0 0.022 0.042 0.062 0.082 0.102 0.122 0.142 0.162 0.182 0.202 0.222 0.242 0.262
3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 1 1 1 1 1 1 1 1 1 1 1 1 1
3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 1 1 1 1 1 1 1 1 1 1 1 1 1
2 2 2 2 2 2 2 2 2 2 2 2 2 2
0.080732 0.079891 0.079334 0.078767 0.078189 0.077600 0.077000 0.076389 0.075767 0.075135 0.074492 0.073838 0.073174 0.072499
242.20 242.09 242.86 243.61 244.34 245.05 245.74 246.42 247.07 247.70 248.31 248.89 249.46 250.00
322.93 327.76 328.85 329.91 330.96 331.98 332.97 333.94 334.88 335.80 336.69 337.55 338.38 339.17
726.59 731.83 734.20 736.52 738.79 741.00 743.16 745.26 747.29 749.27 751.18 753.02 754.80 756.51
403.66 405.54 406.84 408.12 409.37 410.58 411.76 412.91 414.03 415.11 416.16 417.16 418.13 419.07
403.66 403.90 405.18 406.43 407.66 408.85 410.01 411.14 412.23 413.28 414.30 415.29 416.23 417.14
84,265.21 86,214.74 87,495.49 88,803.07 90,138.33 91,502.16 92,895.47 94,319.22 95,774.41 97,262.09 98,783.35 100,339.3 101,931.2 103,560.3
et al., 1993); it can be computed as follows: IS ¼
E½C 1 ð .,.,.,. Þ ; ming E½C g ð .,.,.,. Þ
ð39Þ
where E½C 1 ð .,.,.,. Þ is the minimum total expected cost of generation 1. A larger value of IS is an indication of a better algorithm. Computational time (CT): We measure it by CPU time from the beginning of GA until a termination condition is reached.
As shown in Table 4, we performed the factorial design with four factors, such as population size (P), mutation rate (MR), elitism rate (ER), and selection method, to determine whether they affect the quality of solutions. The details of the factorial design are given in Appendix A. From Table 5, the overall observation is that different population sizes will result in significantly different outputs for the total expected cost, improvement score, convergence score, computational time, and the number of generations. Observations can also be made for an individual response. For a smaller n, the only parameter that affects the total expected cost is the population size. However, for a larger n, there are more factors that affect the total expected cost, and hence these parameters factors should be varied to obtain the optimal solution. For example, for n ¼ 50, the population size, elitism rate, mutation rate, and selection methods should be varied. In addition, it is observed different replications will lead to different computational time requirements. Finally, for a smaller n, the generation numbers of different population sizes are slightly different, but for a larger n, the generation numbers are the same for all factors. We can use four performance indicators (E½C, IS, CS, CT) to compare the GA efficiency. The observations help us analyze the best set of parameters for GA based on the minimum total expected cost in all scenarios, and the best sets of parameters are P ¼30, ER ¼0.4, and MR ¼0.2. From the results of the first set of experiments, it is observed that for a smaller number of items, n, there is no significant difference between the total expected costs obtained from the two selection methods for different mutation rates, elitism rates, and number in population. However, for a larger number of items, n, there is difference in the total expected cost when different mutation rates, elitism rates, and number in population are used, as given in Table 6, where E½CR is the cost determined when RRW is used and E½CT is the cost determined when TS is used. In addition, it is observed that the second method (tournament selection) provides a slightly lower total expected cost. The second set of experiments is implemented for the selected GA parameters given in Table 7, population size, P ¼30, mutation rate, ER ¼0.4, and elitism rate, MR ¼0.2. There are total of
Table 4 Selected values for parameters and selection methods for the first set of experiments. Parameters
Values
Scenarios ðn; KÞ The number of replications Population size (P) Mutation rate (ER) Elitism rate (MR) Selection method Output
(10, 10), (10, 40), (50, 10), (50, 40) 3 10, 20, 30 0.1, 0.2, 0.3, 0.4 0.01, 0.05, 0.1, 0.2, 0.3 RRW, TS E½C, IS, CT, CS, Gen
5 4 2 ¼40 experiments. For both selection methods (RRW and TS), total expected cost, E½C, the improvement score, IS, the CPU time, CT, and the convergence score, CS, are determined. For example, convergence score for RRW and TS are denoted as CSR and CST , respectively, and the results are given in Table 8. With respect to the total expected cost, it is not significantly different in most scenarios. It seems that TS is slightly better than RRW when n increases. As for the improvement score, the larger value the better. It can be seen that, for a smaller value of n, RRW outperforms. However, when n increases, TS outperforms due to a smaller ratio of ISR =IST . In terms of CPU time, when n is small, RRW requires less computational time than that of TS. On the other hand, when n is equal to 50, RRW requires more computational time than that of TS. With respect to the convergence score, the smaller the better. However, observing the ratios of the convergence scores in Table 8, it is inconclusive to determine which method outperforms. In summary, four different measures are used to evaluate the performance of GA. For a smaller number of products, rank roulette wheel performs better based on the total expected cost, improvement score, and CPU time. However, for a larger number of items, tournament selection performs better. Unfortunately, convergence score does not provide any clear conclusion as to which selection method performs better. This indicates that one has to pay attention to the selection method based on the number of items, n, when GA is used to solve the JRP considered in this research. 4.1. Comparison between genetic algorithm and differential evolution algorithm The heuristic approach, differential evolution (DE) algorithm originally presented in Storn and Prince (1997), has several advantages: (a) it can handle different functions such as non-differentiable, nonlinear, and multimodal; (b) it can obtain promising results in a reasonable amount of time by using a vector
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
149
Table 5 Summary of significant p-value of factor to responses of scenarios. Response p-Value of significant factor for each scenario of ðn; KÞ (10, 10)
(10, 40)
(50, 10)
(50, 40)
E½C
Population (0.031)
Population (0.003)
IS CT
Population (0.003) Replication (0.080), Population (0.000), Selection (0.000). Population (0.000) Population (0.047)
Population (0.000) Replication (0.001), Population (0.000), Selection (0.000). Population (0.000) Population (0.002)
Population (0.000), MR (0.052), Selection (0.000). Population (0.000) Replication (0.000), Population (0.000), ER (0.000), Selection (0.000). Population (0.000), Selection (0.000). –
Population (0.000), ER (0.000), Selection (0.000). Population (0.000) Replication (0.000), Population (0.000), MR (0.000), Selection (0.000). Population (0.000), Selection (0.000). –
CS Gen
Table 6 Comparison of the total cost when n¼ 50 and K ¼ $40. ^ER
0.1
0.2
0.3
0.4
0.1
0.2
0.3
0.4
MR
E½CR
0.01 0.05 0.1 0.2 0.3
1,092,778.9 1,093,383.5 1,092,866.0 1,093,100.1 1,093,062.5
1,093,158.4 1,093,229.5 1,092,908.6 1,092,970.1 1,092,834.4
1,093,121.3 1,092,905.4 1,093,300.0 1,092,869.9 1,092,980.8
1,093,051.0 1,092,836.8 1,093,094.8 1,092,896.0 1,092,795.9
1.000044 1.000418 0.999856 1.000307 1.000266
1.000265 1.000513 0.999959 1.000186 0.999889
1.000159 1.000171 1.000410 1.000040 1.000117
1.000154 0.999772 1.000246 1.000185 0.999783
20
0.01 0.05 0.1 0.2 0.3
1,092,734.4 1,092,858.0 1,092,751.6 1,092,702.1 1,092,691.6
1,093,002.3 1,092,733.1 1,092,682.0 1,092,692.3 1,092,757.2
1,092,810.7 1,092,746.3 1,092,900.7 1,092,860.5 1,092,787.2
1,092,807.2 1,092,678.4 1,092,711.9 1,092,698.4 1,092,837.2
0.999932 1.000096 1.000032 0.999950 0.999994
1.000302 1.000067 0.999939 0.999990 1.000047
1.000119 0.999970 1.000199 1.000156 1.000090
1.000088 0.999992 1.000054 1.000015 1.000079
30
0.01 0.05 0.1 0.2 0.3
1,092,898.7 1,092,717.9 1,092,667.9 1,092,665.3 1,092,774.8
1,092,945.3 1,092,742.3 1,092,663.4 1,092,739.8 1,092,700.9
1,093,025.3 1,092,771.8 1,092,659.1 1,092,725.7 1,092,720.9
1,092,864.1 1,092,702.6 1,092,735.1 1,092,684.0 1,092,699.8
1.000200 1.000060 1.000014 1.000012 1.000061
1.000245 1.000070 0.999981 1.000054 1.000030
1.000287 1.000081 1.000004 1.000039 1.000037
1.000176 1.000035 1.000072 0.999972 0.999948
P 10
E½CR =E½CT
population; (c) it has the self-organizing scheme that takes the difference vector of two randomly chosen population vectors to perturb an existing vector; and (d) it has the good convergence properties. The framework of DE is given in Fig. 7. Each step of the framework is described below. 1. Initialization: DE utilizes multi-dimensional parameter vectors to conduct the parallel direct search. Hence, we developed the structure of parameter vectors for two decision variables (mi and bij) as follows: 0
m1
B B m2 B B ⋮ @ mn
b11
b12
⋯
b21
b22
⋯
⋮
⋮
⋯
bn1
bn2
⋯
b1J
Parameters
Values
n K Selection method
5, 10, 20, 30, 50 10, 20, 30, 40 RRW, TS
Eq. (41): ( FðgÞ ¼ 2λ F 0
1
C b2J C C: ⋮ C A bnJ
ð41Þ
λ ¼ e1 Gen þ 1 g ; Gen
Then, we randomly generated the initial values of the vectors according to the uniform distribution U[0, 1] and decoded the values of mi, bij, and Yij as in Section 3.1. 2. Mutation: DE usually adds the scaled difference of two randomly selected parameter vectors to generate a new parameter vector as follows:
νk ðg þ 1Þ ¼ xr1 ðgÞ þF½xr2 ðgÞ xr3 ðgÞ;
Table 7 Selected values for parameters and selection methods for the second set of experiments.
i a r 1 a r 2 ar 3 ;
ð40Þ
where νk ðg þ 1Þ denotes kth new vector for the (g þ1)th generation; and xr1 ðgÞ, xr2 ðgÞ, and xr3 ðgÞ represent r1th, r2th, and r3th target vectors of gth generation, respectively. F is a scale factor that controls the amplification of the differential variation of two target vectors. The decision variables, mi and bij, are mutated based on Eq. (40). The scale factor F is used based on
where F(g) is the scale factor of gth generation, and Gen is the total number of generation. 3. Crossover: In order to increase the diversity of the vectors, crossover generates a trial vector as shown in Eq. (42): uk ðg þ1Þ ¼ ðu1k ðg þ 1Þ; u2k ðg þ 1Þ; …; uzk ðg þ 1ÞÞ; where uzk ðg þ 1Þ ¼
(
νzk ðg þ 1Þ if rand r CR; xzk ðgÞ
if rand 4 CR:
;
z ¼ 1; 2; …; Z;
ð42Þ
ð43Þ
where rand is a uniform random number within the interval [0,1], CR is the user-defined crossover rate, Z denotes the number of dimension of vectors. The decision variables, mi and bij, are crossovered according to Eqs. (42) and (43). 4. Selection: In this step, the trial vector will be compared with the target vector based on the total expected cost function f. The
150
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
Table 8 Comparison of performance measures and for selected GA parameters. n
E½CR
ISR
CTR
CSR
GenR
E½CR =E½CT
ISR =IST
CT R =CT T
CSR =CST
GenR =GenT
5 10 20 30 50
88,035.3 206,909.5 446,537.9 645,911.6 1,092,277.0
1 1.000114 1.000134 1.000875 1.003289
0.5332 0.111575 1.034626 2.181036 8.015088
1 1.000022 1.000011 1.000098 1.000720
10 11 61 81 101
1 1 1 1 1.00004
1 0.999967 0.999934 1.000118 0.997390
4.299147 0.723527 0.669337 0.728164 1.048020
1 1.000006 0.999996 1.000033 0.999893
1 1.1 1 1 1
20
5 10 20 30 50
88,159.1 207,027.0 446,656.8 646,061.5 1,092,402.0
1 1.000079 1.000233 1.000648 1.006614
0.082572 0.106939 1.056045 2.196840 8.045402
1 1.000017 1.000014 1.000092 1.001205
10 10 61 81 101
1 1 1 1.000014 1.000028
1 1.000079 0.999935 1.000015 1.002205
0.668619 0.677760 0.687305 0.735720 1.049723
1 1.000017 0.999999 1.000045 1.000536
1 1 1 1 1
30
5 10 20 30 50
88,278.01 207,142.0 446,774.5 646,182.6 1,092,573.0
1 1.000085 1.00027 1.001048 1.005727
0.082319 0.120387 1.052576 2.194502 8.072410
1 1.000018 1.000015 1.000111 1.000898
10 10 61 81 101
1 1 1 1 1.000038
1 1 1.000091 1.000025 0.999572
0.661697 0.759408 0.684184 0.733724 1.053705
1 1.000004 0.999994 1.000051 0.999889
1 1 1 1 1
40
5 10 20 30 50
88,375.6 207,254.8 446,891.0 646,301.3 1,092,653.0
1 1.000002 1.000152 1.000578 1.006092
0.084821 0.119709 1.034572 2.199863 8.007101
1 1 1.000025 1.000078 1.000812
10 10 61 81 101
1 0.999999 1 0.999956 0.999968
1 0.999953 0.999971 0.999972 1.001295
0.685297 0.748792 0.667291 0.730768 1.056270
1 0.999988 1.000013 1.000025 1.000240
1 1 1 1 1
K 10
be seen that, for a larger value of n, GA obtained a lower total expected cost but required more computational time than that of DE. Nonetheless, overall, the total expected cost obtained from GA is 2.31% less than that of DE. Similarly, computational time (CT) required by GA is 0.52% less than that of DE. Hence, GA is better in terms of these two measures. IS and CS scores determine how well the solutions are improved in subsequent generations. Table 9 shows that most IS values obtained from DE are larger than those obtained from GA, and hence DE is better. Similarly, CS scores obtained from DE are close to one and it implies that the average value of fitness functions is not different from the best value. In addition, most CS scores obtained from DE are slightly larger than those obtained from GA, and hence DE is better. After all, the main objective of the comparison is the total expected cost, and hence the GA is preferred.
5. Conclusion Fig. 7. DE algorithm.
vector with the smaller value of the total expected cost will be put into the next generation ( uk ðg þ 1Þ; f ðuk ðg þ1ÞÞ r f ðxk ðgÞÞ; xk ðg þ1Þ ¼ ð44Þ Otherwise: xk ðgÞ The decision variables, mi and bij, used Eq. (44) to select the child vectors of the next generation. 5. Stopping criterion: We used the same stopping criterion as in the GA. The proposed JRP problem was solved using the DE algorithm and then the results obtained from DE and GA with selection method of rank roulette wheel are compared as shown in Table 9. Experiments are implemented for different values of K and n. The recorded outputs are the total expected cost, improvement score, computational time, and the number of generations. In order to compare the results, differences of outputs from both GA and DE are presented in the last five columns of Table 9. For example, the difference in the total expected cost for K ¼ 10 and n ¼5 is calculated as fð88; 035:3–90; 247:18Þ=88; 035:3g 100 ¼ 2:51%. It can
Resource constraints are always inherent with a JRP; there have been several studies incorporating a few of them. This research simultaneously takes into account not only the constraints often addressed in the literature, such as budget constraint and capacity constraint, but also dependency of items and defective items in each shipment. This adds a higher dimension to the problem and captures the nature of a real JRP problem existing in many retail companies. In order to investigate the ordering policy, the optimal family cycle length, the best set of order frequencies, the order quantity, and the total expected cost per unit time are determined. A two-dimensional genetic algorithm (GA) is proposed to solve the problem. An empirical study is conducted to determine the best GA parameters and selection method. The results also show that a higher percentage of defective items leads to a lower optimal family cycle length, a higher order quantity for each item, and a higher total expected cost per unit time. The JRP is also solved by using the differential evolution (DE) algorithm and results obtained from DE are compared with those obtained from GA. Based on the minimum total expected cost of the JRP, GA is preferred. This research can be extended in several ways. There are a number of ways to add more aspects to the investigated problem.
Table 9 The comparison between GA and DE. n
GA-RRW
DE
Difference (%)
E½C
IS
CT
CS
Gen
E½C
IS
CT
CS
Gen
E½C
IS
CT
CS
Gen
5 10 20 30 50
88,035.3 206,909.5 446,537.9 645,911.6 1,092,277.0
1 1.000114 1.000134 1.000875 1.003289
0.5332 0.111575 1.034626 2.181036 8.015088
1 1.000022 1.000011 1.000098 1.00072
10 11 61 81 101
90,247.18 211,120.9 454,091.3 654,215.5 1,134,355
1 1.01712 1.015416 1.019667 1.039289
0.124276 0.177396 1.553928 2.806978 2.030138
1 1.007874 1.005948 1.005541 1.010408
10 11 61 81 101
2.51 2.04 1.69 1.29 3.85
0.00 1.70 1.53 1.88 3.59
76.69 58.99 50.19 28.70 74.67
0.00 0.79 0.59 0.54 0.97
0 0 0 0 0
20
5 10 20 30 50
88,159.1 207,027.0 446,656.8 646,061.5 1,092,402.0
1 1.000079 1.000233 1.000648 1.006614
0.082572 0.106939 1.056045 2.19684 8.045402
1 1.000017 1.000014 1.000092 1.001205
10 10 61 81 101
90,298.52 211,171.6 453,001.9 654,282.5 1,159,125
1 1.011971 1.011447 1.019677 1.024488
0.121747 0.161616 1.554122 2.828463 2.089074
1 1.005273 1.002913 1.006586 1.005818
10 10 61 81 101
2.43 2.00 1.42 1.27 6.11
0.00 1.19 1.12 1.90 1.78
47.44 51.13 47.16 28.75 74.03
0.00 0.53 0.29 0.65 0.46
0 0 0 0 0
30
5 10 20 30 50
88,278.0 207,142.0 446,774.5 646,182.6 1,092,573.0
1 1.000085 1.00027 1.001048 1.005727
0.082319 0.120387 1.052576 2.194502 8.07241
1 1.000018 1.000015 1.000111 1.000898
10 10 61 81 101
90,349.37 210,157.8 454,202.4 655,691.3 1,139,811
1 1.012068 1.014232 1.032538 1.040375
0.122499 0.163604 1.552599 2.854533 2.031563
1 1.003739 1.002678 1.010285 1.008621
10 10 61 81 101
2.35 1.46 1.66 1.47 4.32
0.00 1.20 1.40 3.15 3.45
48.81 35.90 47.50 30.08 74.83
0.00 0.37 0.27 1.02 0.77
0 0 0 0 0
40
5 10 20 30 50
88,375.7 207,254.8 446,891.0 646,301.3 1,092,653.0
1 1.000002 1.000152 1.000578 1.006092
0.084821 0.119709 1.034572 2.199863 8.007101
1 1 1.000025 1.000078 1.000812
10 10 61 81 101
90,399.72 210,214.1 453,135.5 6,557,58.4 1,133,157
1 1.015391 1.016618 1.035515 1.041107
0.123717 0.159229 1.554821 2.877044 2.031391
1 1.006917 1.006004 1.01086 1.008839 Average
10 10 61 81 101 2.31
2.29 1.43 1.40 1.46 3.71 1.70
0.00 1.54 1.65 3.49 3.48 12.99
45.86 33.01 50.29 30.78 74.63 0.52
0.00 0.69 0.60 1.08 0.80 0
0 0 0 0 0
10
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
K
151
152
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
For example, incorporating quantity discounts in this particular model could be an extension. Multiple suppliers could also be considered in future research. Moreover, instead of considering complete screening of the shipment, an acceptance sampling approach can be implemented. In terms of solution methodology, differential evolution algorithm can be further improved and then can be compared with the proposed GA. In this research, we considered that the percentage of defective of an item follows a uniform distribution. However, the percentage of defective of an item can follow any other distribution and such consideration can be investigated in the future research.
Acknowledgment The authors would like to thank Natural Sciences and Engineering Research Council of Canada and National Nature Science Foundation of China for funding this research under the Grant number of 71402149.
Appendix A. The details of the factorial design Each factor is varied to the levels that they will more likely affect the results. For example, population size is varied at three levels: 10, 20 and 30. Elitism rate is varied at four levels: 0.1, 0.2, 0.3, and 0.4. Mutation rate is varied at five levels: 0.01, 0.05, 0.1, 0.2, and 0.3, and it is noted that if elitism rate is more than 0.3, the result converges very quickly, and two selection methods are considered: Rank Roulette Wheel (RRW) and Tournament Selection (TS). Then, GA is repeated three times for four combinations of items and the ordering cost, i.e., (n,K) ¼{(10,10), (10,40), (50,10), (50,40)}. The analysis of variance is performed for five responses: the total expected cost ðE½CÞ, improvement score (IS), cpu time (CT), convergence score (CS), and the number of generations (Gen). Hence, for four combinations of ðn; KÞ, there will be 3 3 4 5 2 4 ¼ 1440 scenarios to be solved, but each scenario has five responses and hence the total number of responses is 7200. To be efficient, we performed the main effect analysis and determined which factor affects the response and the results are concluded in Table 5. The main effect analysis is tested with 95% confidence interval (i.e., α ¼0.05) and p-value is given within the parentheses. The p-value indicates if the outputs are significantly different for different levels of the factors.
References Anderson, C., Jones, K.F., Ryan, J., 1991. A two-dimensional genetic algorithm for the ising problem. Complex Syst. 5 (3), 327–333. Beasley, D., Bull, D.R., Martin, R.R., 1993. An overview of genetic algorithms: Part 1, fundamentals. Univ. Comput. 15 (2), 58–69. Bortfelldt, A., 2006. A genetic algorithm for the two-dimensional strip packing problem with rectangular pieces. Eur. J. Oper. Res. 172 (3), 814–837. Bortfelldt, A., Winter, T., 2009. A genetic algorithm for the two-dimensional knapsack problem with rectangular pieces. Int. Trans. Oper. Res. 16 (6), 685–713. Cárdenas-Barrón, L.E., 2000. Observation on: “economic production quantity model for items with imperfect quality” [Int. J. production economics 64 (2000) 59–64]. Int. J. Prod. Econ. 67 (2), 201.
Cha, B.C., Moon, I.K., 2005. The joint replenishment problem with quantity discounts under constant demand. OR Spectr. 27, 569–581. De Jong, K., 1980. Adaptive system design: a genetic approach. IEEE Trans. Syst. Man Cybern. 10 (9), 566–574. Goyal, S.K., 1974. Determination of optimum packaging frequency of items jointly replenished. Manag. Sci. 21 (4), 436–443. Goyal, S.K., 1975. Analysis of joint replenishment inventory systems with resource restriction. Oper. Res. Q. 26 (1), 197–203. Goyal, S.K., Cárdenas-Barrón, L.E., 2002. Note on: economic production quantity model for items with imperfect quality—a practical approach. Int. J. Prod. Econ. 77, 85–87. Hariga, M., 1994. Two new heuristic procedures for the joint replenishment problem. J. Oper. Res. Soc. 45 (4), 463–471. Herbert, E.A., Dowsland, K.A., 1996. A family of genetic algorithms for the pallet loading problem. Ann. Oper. Res. 63 (3), 415–436. Hong, S.P., Kim, Y.H., 2009. A genetic algorithm for joint replenishment based on the exact inventory cost. Comput. Oper. Res. 36, 167–175. Jia, W., Jiang, L., Li, H., Zheng, G., Li, X., 2009. Genetic algorithm for band gap optimization under light line in two-dimensional photonic crystal slab. Opt. Appl. 39 (3), 481–488. Kaspi, M., Rosenblatt, M.J., 1983. An improvement of Silver's algorithm for joint replenishment problem. IIE Trans. 15 (3), 264–267. Kaspi, M., Rosenblatt, M.J., 1991. On the economic ordering quantity for jointly replenished items. Int. J. Prod. Res. 29, 107–114. Khouja, M., Goyal, S.K., 2008. A review of the joint replenishment problem literature: 1989–2005. Eur. J. Oper. Res. 186 (1), 1–16. Khouja, M., Michalewicz, Z., Satoskar, S.S., 2000. A comparison between genetic algorithm and the RAND method for solving the joint replenishment problem. Prod. Plan. Control 11 (6), 556–564. Kruse, R., Borgelt, C., Klawonn, F., Moewes, C., Steinbrecher, M., Held, P., 2013. Computational Intelligence: A Methodological Introduction. Springer-Verlag, London. Maddah, B., Jaber, M.Y., 2008. Economic order quantity for items with imperfect quality: revisited. Int. J. Prod. Econ. 112, 808–815. Moon, I.K., Cha, B.C., 2006. The joint replenishment problem with resource restriction. Eur. J. Oper. Res. 173 (1), 190–198. Olsen, A.L., 2008. Inventory replenishment with interdependent ordering costs: an evolutionary algorithm solution. Int. J. Prod. Econ. 113, 359–369. Ono, T., Ikeda, T., 1998. Optimization of two-dimensional guillotine cutting by genetic algorithms. In: The 6th Proceedings of The European Congress on Intelligent Techniques and Soft Computing, (EUFIT 1998), 7–10 September 1998, Aachen, Germany, pp. 450–454. Papachristos, S., Konstantaras, L., 2006. Economic ordering quantity models for items with imperfect quality. Int. J. Prod. Econ. 100, 148–154. Porteus, E., 1986. Optimal lot sizing, process quality improvement, and setup cost reduction. Oper. Res. 34 (1), 137–144. Robinson, P., Narayanan, A., Sahin, F., 2009. Coordinated deterministic dynamic demand lot-sizing problem: a review of models and algorithms. Omega 37 (1), 3–15. Rosenblatt, M., Lee, H., 1986. Economic production cycles with imperfect production processes. IIE Trans. 18 (1), 48–55. Salameh, M.K., Jaber, M.Y., 2000. Economic production quantity model for items with imperfect quality. Int. J. Prod. Econ. 64, 59–64. Shu, F., 1971. Economic ordering frequency for two items jointly replenished. Manag. Sci. 17, 406–410. Silver, E.A., 1975. Modifying the economic order quantity (EOQ) to handle coordinated replenishment of two or more items. Prod. Inventory Manag. 16, 26–38. Storn, R., Prince, K., 1997. Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. van Eijis, M.J.G., Heuts, R.M.J., Kleijnen, J.P.C., 1992. Analysis and comparison of two strategies for multi-item inventory systems with joint replenishment costs. Eur. J. Oper. Res. 59 (3), 405–412. Wahab, M.I.M., Jaber, M.Y., 2010. Economic order quantity model for items with imperfect quality, different holding costs, and learning effects: a note. Comput. Ind. Eng. 58 (1), 186–190. Wahab, M.I.M., Mamun, S., Ongkunaruk, P., 2011. EOQ models for a coordinated two-level international supply chain considering imperfect items and environmental impact. Int. J. Prod. Econ. 134, 151–158. Wang, L., He, J., Wu, D., Zeng, Y.-R., 2012. A novel differential evolutional algorithm for joint replenishment problem under interdependence and its application. Int. J. Prod. Econ. 135, 190–198. Wildeman, R.E., Frenk, J.B.G., Dekker, R., 1997. An efficient optimal solution method for the joint replenishment problem. Eur. J. Oper. Res. 99 (2), 433–444.
Contents lists available at ScienceDirect
Int. J. Production Economics journal homepage: www.elsevier.com/locate/ijpe
A genetic algorithm for a joint replenishment problem with resource and shipment constraints and defective items P. Ongkunaruk a, M.I.M. Wahab b,n, Y. Chen c a
Department of Agro-Industrial Technology, Kasetsart University, 50 Ngam wong wan Road, Lad Yao, Chatuchak, Bangkok 10900, Thailand Department of Mechanical and Industrial Engineering, Ryerson University, 350 Victoria Street, Toronto, ON, Canada 5MB 2K3 c School of Economics and Management, Southwest Jiaotong University, Chengdu, Sichuan 610031, PR China b
art ic l e i nf o
a b s t r a c t
Article history: Received 3 June 2015 Accepted 18 February 2016 Available online 2 March 2016
A joint replenishment problem (JRP) is presented to determine the optimal reordering policy for multiple items with defective quantity and several restrictions such as shipment constraint, budget constraint, and transportation capacity constraint. The objective is to minimize the total expected cost per unit time. A two-dimensional genetic algorithm (GA) is provided to determine an optimal family cycle length and the reorder frequencies. A numerical example is presented and the results are discussed including the effect of defective items on the ordering policy. Extensive computational experiments are performed to test the performance of the GA. The JRP was also solved by a differential evolution (DE) algorithm and the results obtained from DE were compared with those obtained from GA. Crown Copyright & 2016 Published by Elsevier B.V. All rights reserved.
Keywords: Joint replenishment Defective items Shipment constraint Genetic algorithm Differential evolution algorithm
1. Introduction A joint replenishment problem (JRP) is defined as the coordination of the replenishment of a group of items that may be ordered jointly from a single supplier (Goyal, 1974). One of the advantages of jointly placing order for a group of items is sharing the fixed cost associated with order, and hence the total cost is reduced when compared to the cost of individually placing order for each item (Moon and Cha, 2006). Moreover, a JRP is a good practice to reduce transportation cost, when a vehicle is used to transport multiple items simultaneously. Similarly, the supplier's set-up cost can be reduced by manufacturing multiple items using the same facility. One of the disadvantages of the JRP is increasing the average inventory level and system control costs. In order to solve the basic JRP, family cycle time and the frequency of replenishing individual items have to be determined while minimizing the total cost, which is the sum of ordering and holding costs. Most of the studies related to the basic JRP have usual assumptions such as deterministic and constant demand, linear holding cost, no shortage, and no quantity discount (e.g., Khouja and Goyal, 2008; Robinson et al., 2009). The JRP has received considerable attention over the last few decades and different solution approaches have been proposed to solve the JRP. For example, Shu (1971), Goyal (1974), and Silver n
Corresponding author: Tel.: þ 1 416 979 5000x2670; fax: þ1 416 979 5265. E-mail addresses: [email protected] (P. Ongkunaruk), [email protected] (M.I.M. Wahab), [email protected] (Y. Chen). http://dx.doi.org/10.1016/j.ijpe.2016.02.012 0925-5273/Crown Copyright & 2016 Published by Elsevier B.V. All rights reserved.
(1975) proposed a simple procedure to determine the order quantity for multiple items in a JRP. Moreover, many heuristics have been developed to solve different versions of JRP, namely, basic JRP, JRP with quantity discount, and constrained JRP using various methods including direct and indirect grouping strategies, genetic algorithm (GA), evolutionary algorithm (EA), RAND, QD– RAND, and C-RAND (e.g., Kaspi and Rosenblatt, 1983, 1991; van Eijis et al., 1992; Hariga, 1994; Wildeman et al., 1997; Khouja et al., 2000; Cha and Moon, 2005; Moon and Cha, 2006; Olsen, 2008; Hong and Kim, 2009). However, among many heuristics, GA and EA have been proved to be effective for solving the JRP (Olsen, 2008). Moreover, it was concluded that solutions obtained from GA combined with a local search algorithm and RAND are very close to the global optimum (Hong and Kim, 2009). Recently, Wang et al. (2012) proposed a differential evolution (DE) algorithm based on Storn and Prince (1997) to solve a JRP and claimed that the DE algorithm can solve a JRP more effectively than the EA used in Olsen (2008). Hence, in this research, we used GA and DE, which are the best algorithms existing in the literature, to solve a JRP. With respect to GA literature, we used a two dimensional chromosome representation. There are several studies on twodimensional GA for solving packing problem (e.g., Herbert and Dowsland, 1996; Bortfelldt, 2006), cutting problem (e.g., Ono and Ikeda, 2008), and knapsack problem (e.g., Bortfelldt and Winter, 2009). In addition, the two-dimensional GA was also implemented in different areas of research, for example, ising problem (Anderson et al., 1991) and photonic crystal slab (Jia et al., 2009). With respect to DE, Storn and Prince (1997) proposed a novel DE to
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
solve a complex continuous non-linear functions and we used a similar algorithm to solve the JRP. Most of the studies on JRP assumed that items included in an order can be shipped together without any restrictions. However, in practice, there are restrictions or dependency among the items. As a result, though an order is placed for certain items at the same time, they cannot be shipped together in the same truck. For example, chemical or fragrant products should not be shipped together with dry food or food ingredients, since the odor will be absorbed by the dry food or food ingredients. Another example is ambient products, which are not shipped with fresh produce and ready to eat meals due to different storage temperatures. Recently, several studies (e.g., Olsen, 2008; Wang et al., 2012) considered such interdependency of the items in a JRP. However, these studies assumed that such interdependency increases minor ordering costs to maintain different requirements of the prohibited items in the same truck and then the JRP was solved by using direct and indirect groupings. Let us consider a modern retail business company, which operates in Thailand, having restrictions in its shipment policy to maintain quality and food safety. This company does not allow transporting prohibited items in the same truck at all. This retail company plans its ordering policy to minimize its total inventory related cost considering quality and safety requirements of the company. Then, the company provides the ordering policy to a third-party logistics provider and it ships the items for a fixed transportation cost. Hence, the JRP models in Olsen (2008) and Wang et al. (2012) cannot be applied to this retailer. In order to solve the retailer's JRP, in this paper, without changing minor setup costs of items having restrictions, shipment restrictions or interdependency of the items are included in the constraints and then a constrained JRP is solved. In addition, there are several other constraints besides shipping restrictions or interdependency of items. For example, Goyal (1975) introduced a JRP with budget constraint. Khouja et al. (2000) addressed the issue of maximum load constraint in a JRP. Later, Moon and Cha (2006) also studied a JRP with a budget constraint. However, the existing literature on JRP does not consider such constraints simultaneously in a JRP. In this research, we included budget constraint and transportation capacity constraint with interdependency of the items and then solved a constrained JRP. Another important aspect in a JRP is defective items. Wahab et al. (2011) pointed out that many Thai small and medium enterprises implemented total quality management, however buyers still sometimes find defective items in the shipments, and it is also possible that some items may also get damaged during shipment. There are several studies dealing with defective items in a single item replenishment such as EOQ, for example, Porteus (1986), Rosenblatt and Lee (1986), Salameh and Jaber (2000), Cárdenas-Barrón (2000), Goyal and Cárdenas-Barrón (2002), Papachristos and Konstantaras (2006), Maddah and Jaber (2008), and Wahab and Jaber (2010). The defective items in shipments will affect the ordering policy. For example, if there are more defective items in the shipment, the order quantity should be higher to compensate for those defective items. Hence, this study also considered the percentage of defective items in the JRP to investigate its effect in the ordering policy. The purpose of this research is to model a constrained JRP for multiple items, where each item has a certain percentage of defective item in each shipment, with constraints such as dependency among items (i.e., some items have restrictions that they cannot be shipped together in the same truck), budget constraint, and transportation capacity constraint. Such a JRP exists in many real business applications. For example, this problem can be found in a retailer that delivers products from its distribution center to retail stores. The JRP is solved by using a two-dimensional GA and DE. The results obtained from both GA and DE are compared.
143
This paper is organized as follows: Section 2 presents the model, notations, and decision variables. Section 3 presents the solution methodology, a two-dimensional GA, and its steps. Section 4 illustrates numerical examples and also presents the comparison between GA and DE. Section 5 concludes the paper.
2. The model This section first presents indices, parameters, and decision variables of the model and then describes the JRP step by step to develop the model. 2.1. Indices, parameters, and decision variables Indices i an index of an item, where i ¼ 1; 2; …; n j an index of a truck that is used to ship items, where j ¼ 1; 2; …; J l an index of a set representing items that can not be shipped together, where l ¼ 1; 2; …; L Parameters n the number of items J the number of trucks L the number of sets of items that cannot be shipped together Δj the capacity of truck j Di Qi ai
the the the the
demand rate for item i order quantity for item i handling cost for item i unit screening cost for item i
μi νi the lost per defective unit for item i pi the percentage of defective of item i xi the screening rate for item i ti the screening time for item i hi the holding cost per unit of item i per unit time ci the unit variable cost for item i K the ordering cost per order B the amount of capital that can be invested Sl the lth set of items that cannot be shipped together Decision variables T time interval between replenishment for a set of items mi the integer number of intervals that the replenishment takes place for item i bij the proportion of item i transported by truck j ( 1 if item i is ordered and delivered by truck j; Yij 0; otherwise:
2.2. The joint replenishment problem Item i has demand Di per year. It is assumed that demand for each item is deterministic and constant, and shortages are not allowed. Item i is ordered in every mith replenishment cycles, where mi is a positive integer, meaning the replenishment quantity of item i will last for miT period of time. This also indicates that the cycle time for item i is miT. An important aspect that has been addressed in this problem is the restriction in shipping. That is, because of the nature of the products and requirements, some items cannot be shipped together. For example, chemical or fragrant products should not be shipped in the same truck as dry food or raw material. Though
144
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
these products were ordered at the same time, they are shipped in different trucks with other products that do not have any conflict. In other words, several items can be shipped simultaneously by a truck if there is no restrictions in shipping them together. In this scenario, depending on capacity requirements for shipping, sometimes an ordered quantity of an item may not be shipped by a single truck and it can be divided and shipped by more than one truck at the same time. A joint ordering cost K is incurred every T period of time. For example, the cost of placing an order for multiple items via fax, email, or phone, and this cost is independent of types of items and their order quantities. The handling cost ai incurs after receiving item i for every miT period of time. Therefore, the total setup cost per unit time is Kþ
Pn
i¼1
PJ j¼1
ai bij Y ij mi
T
;
n X
c i Di : 1 E½p i i¼1
Q i ð1 E½pi Þ Z Di t i
ð9Þ
The total lost of defective units, due to markdown, is given by n X νi Di E½p i
i¼1
1 E½pi
:
ð10Þ
The total expected cost per unit time for the JRP is the summation of all aforementioned costs and it is given as follows: 0
E½CðT; m0i s; bij s; Y 0ij sÞ ¼
n X Di ðci þ μi þ νi E½pi Þ 1 E½pi i¼1
Di : xi
ð2Þ
ð3Þ
ð5Þ
ð6Þ
By substituting for Qi and ti, the above equation can be expressed as: Di mi T D2 m TE½pi þ i i : 2 xi ð1 E½pi Þ2
þ
PJ j¼1
ai bij Y ij mi
T " # D i mi T D2 m TE½pi þ i i þ hi : 2 xi ð1 E½pi Þ2 i¼1
ð11Þ
This cost function needs to be minimized subject to capacity, budget, and shipping constraints as follows: n X Di mi bij Y ij T i¼1
ð1 E½pi Þ
r Δj ;
8j
ð12Þ
n X D i mi c i T r B; ð1 E½pi Þ i¼1
ð13Þ
bij r Y ij ;
ð14Þ
8 i; j
0 r bij r 1;
8 i; j
ð15Þ
bij ¼ 1;
8i
ð16Þ
Y ij r 1;
8 j; l
ð17Þ
J X
X i A Sl
T Z0
ð18Þ
mi A f1; 2; …g; Y ij A f0; 1g;
ð4Þ
where Qi consists of good and defective items, and their average inventory level is given by Q i ð1 E½pi Þ Q i E½pi t i þ 2 mi T
i¼1
j¼1
Since only good items can be used to meet the demand during miT periods of time, the replenishment quantity Qi for item i is given by Q i ð1 E½pi Þ ¼ Di mi T;
Pn
n X
Substituting t i ¼ Q i =xi , it becomes
Ii ¼
μi Di : 1 E½pi i¼1
ð1Þ
Moreover, the following condition must be satisfied for no-shortage:
Ii ¼
n X
Kþ
in this expression, the second term indicates that the handling cost is charged for each fraction bij of the order quantity shipped in each vehicle. In most industrial applications, due to the damage in transit and/or imperfect production process, some defective items present in a shipment. It is assumed that pi is the percentage of defective for item i in each replenishment, where pi is a random variable. Once the replenishment is received, it is screened at the rate of xi, which is greater than Di, during the time t i ¼ Q i =xi . When the defective items are segregated by the screening process, those defective items are immediately sold for a markdown price, and hence defective units of item i are retired from the system at time ti. For item i, the expected number of good units in each replenishment is Q i ð1 E½pi Þ, and only good units can be used to meet the demand in each cycle. Hence, the total purchasing cost per unit time is
E½pi r1
Since μi is the screening cost per unit, the total screening cost per unit time is given by
ð7Þ
A holding cost of hi for item i is incurred per unit time. Therefore, the total expected holding cost per unit time is " # n X DmT D2 m TE½pi hi i i þ i i : ð8Þ 2 xi ð1 E½pi Þ2 i¼1
8i
ð19Þ
8 i; j
ð20Þ
Eq. (12) is the capacity constraint. The total shipment size in a truck cannot be more than its capacity. Eq. (13) is the budget constraint. The total purchasing per order should not be more than the budget. Inequality (14) is to ensure that if Yij is zero, then bij is zero. Inequalities in (15) are the boundary constraints of bij. Eq. (16) implies that item i can be fully or partially delivered by truck j. Inequality (17) is the shipment requirement constraint. Items in set Sl cannot be shipped simultaneously in the same truck. Eq. (18) implies that the time interval must be positive. Eq. (19) implies that mi must be positive integer that is greater than or equal to one. Eq. (20) implies that Yij is binary. The objective is to determine the decision variables, mi, T, bij, and Yij by minimizing the total expected cost per unit time given in Eq. (11). Therefore, taking the first derivative of the total expected cost with respect to family cycle length T and setting it to zero, we obtain Kþ
0
dE½CðT; m0i s; bij s; Y 0ij sÞ dT
Pn
¼
i¼1
PJ j¼1
ai bij Y ij mi
T2 # Di mi D2i mi E½pi þ þ hi ¼ 0: 2 xi ð1 E½pi Þ2 i¼1 n X
"
ð21Þ
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
145
After rearranging the above equation, we obtain the following solution for T and it is denoted as T0. vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P P ai bij Y ij u K þ ni¼ 1 Jj ¼ 1 u u mi " #: T0 ¼ u ð22Þ 2 uP t n h Di mi þ Di mi E½pi i¼1 i 2 xi ð1 E½pi Þ2 From the capacity constraint, Eq. (12), it can be written as: T 1 ¼ min
j ¼ 1;…;J Pn i¼1
Δj Di mi bij Y ij ð1 E½pi Þ
:
ð23Þ
From the budget constraint, Eq. (13), it can also be written as: T2 ¼
B : Di mi ci i¼1 ð1 E½pi Þ
ð24Þ
Pn
0
Finally, for a given set of m0i s, bij s, and Y 0ij s, the optimal cycle time Tn can be determined as follows: T n ¼ minfT 0 ; T 1 ; T 2 g:
ð25Þ 0 bij s,
Y 0ij s
and that The next step is to determine the best set of m0i s, minimize the total expected cost per unit time given in Eq. (11). In 0 order to determine m0i s, bij s, and Y 0ij s, GA is employed in the next section.
3. Genetic algorithm Genetic algorithm starts with an initial set of random solutions called a population. Each solution in the population is called a chromosome. Each chromosome composes a set of genes. All genes are generated by the random number generator of uniformly distributed random number between 0 and 1. For each generation, the chromosomes are evaluated using some measures of fitness. The algorithm is presented in Fig. 1 and explained as follows: 3.1. GA chromosome representation Fig. 1. GA.
The basic cycle time T, n integer values of m0i s, n J values of and n J values of Y 0ij s have to be decided for solving the JRP. In order to determine the values of those decision variables, the following steps are given with examples. 0 bij s,
1. Encode the chromosome that has n ðJ þ 1Þ genes. Each gene represents a uniform random variable between 0 and 1. An example is given in Fig. 2 for i¼ 5 and j¼ 2. The first column represents the values of m0i s and the rest represent the values of 0 bij s. 2. Decode the chromosome that represents the solution of the JRP as follows: GeneðiÞð1Þ that represents mi can be decoded in the following way: UB LB mi ¼ mLB i þ⌊ðmi mi þ1Þ GeneðiÞð1Þc UB
8 i ¼ 1; 2; …; n;
ð26Þ
LB
where the upper bound, mi , and a the lower bound, mi , of mi can be determined from Eqs. (33) and (34), respectively. Then, there are two cases to be considered as follows: (a) If there is no restriction that items can be shipped in the same 0 truck, GeneðiÞð2Þ with n ðJ þ 1Þ genes that represent bij s can be decoded in the following way: bi1 ¼ GeneðiÞð2Þ
8 i ¼ 1; 2; …; n:
ð27Þ
Fig. 2. Encoding chromosome with n ðJ þ 1Þ genes.
bij ¼
1
j1 X
! bik
k¼1
GeneðiÞðj þ1Þ
biJ ¼ 1
J1 X
bij
8 i ¼ 1; 2; …; n; j ¼ 2; …; J 1:
8 i ¼ 1; 2; …; n:
ð28Þ
ð29Þ
j¼1
Next, determine Y 0ij s if bij 4 0; then Y ij ¼ 1; otherwise, Y ij ¼ 0. An example is given in Fig. 3.
146
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
(b) If there is restriction that some items cannot be shipped together in the same truck, then Y 0ij s, which do not violate Eq. (17), are randomly determined. Next, Gene ðiÞð2Þ with n 0 ðJ þ1Þ genes that represent bij s can be decoded in the following way: bi1 ¼ GeneðiÞð2Þ Y i1
bij ¼
1
j1 X
8 i ¼ 1; 2; …; n:
GeneðiÞðjþ 1Þ
k¼1
Y ij
biJ ¼ 1
8 i ¼ 1; 2; …; n; j ¼ 2; …; J 1:
J1 X
bij
8 i ¼ 1; 2; …; n:
ð31Þ
ð32Þ
j¼1
An example is given in Fig. 4. It is well-known that if the cost is minimum at mi, then it is lower than the cost at mi 1 or at mi þ 1 . Hence, in order to minimize the cost function given in Eq. (11), the lower bound of mi, miLB, and the upper bound of mi, miUB, can be defined from the optimality condition as follows (e.g., Goyal, 1974): "
UB mUB i ðmi 1Þ r
hi T 2min
ai D2i E½pi
Di þ 2 xi ð1 E½pi Þ2
UB # r mUB i ðmi þ 1Þ;
8 i;
ð33Þ "
LB mLB i ðmi 1Þ r
hi T 2max
ai Di D2i E½pi þ 2 xi ð1 E½pi Þ2
i
ð36Þ
i
ð30Þ
! bik
9 8 > > > > > > >vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi> > > = > i ¼ 1;…;n>u > Di E½pi > > t Di > > > ; : hi 2 þ x ð1 E½p Þ2 >
LB # r mLB i ðmi þ 1Þ;
8 i;
ð34Þ where Tmax can be estimated from Eq. (22) by substituting mi equal to 1; Tmin can be obtained from the standard EOQ for an individual item. They can be expressed as given below: vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u P P u K þ ni¼ 1 Jj ¼ 1 ai bij " #: ð35Þ T max ¼ u u uPn Di D2i E½pi t þ i ¼ 1 hi 2 xi ð1 E½pi Þ2
3. Determine the optimal basic cycle time, Tn, from Eq. (25) by substituting a given set of m1 ; m2 ; …; mn , b11 ; …; bnJ , and Y 11 ; …; Y nJ . 4. Check the feasibility of the solutions by substituting them in Eqs. (12), (13), and (17). 5. Evaluate the fitness function. Calculate the total expected cost per unit time by using the objective function given in Eq. (11). Then, the fitness function is the summation of the objective function and penalty cost. There are three constraints: vehicle capacity, budget, and dependency of items. If the chromosome violates any constraint, then a larger penalty cost is added to the fitness function. Hence, the probability of being selected to be a parent is reduced. Otherwise, no penalty cost is added to the fitness function. 3.2. Selection We implement GA to select the best method by applying two selection methods: (i) a combination of rank and roulette wheel selection (Kruse et al., 2013), and (ii) an unbiased tournament. (i) Using the first method, initially the fitness function is ranked in the ascending order, and then the probability of selection is calculated. The probability of selection is determined by dividing its rank by the total summation of the rank. (ii) Unbiased tournament randomly divides the chromosomes into groups and selects the best solution from each group. 3.3. GA operators GA operators create a new generation from parent chromosomes. There are three GA generators as follows: (i) Elitism: In order to maintain the best solutions from the last generation to the next generation, a certain number of best chromosomes will be copied to the next generation with a specific probability (Beasley et al., 1993; Kruse et al., 2013). (ii) Crossover: The crossover of the selected parents will produce the child for the next generation. We apply horizontal crossover by randomly selecting any row of two parent chromosomes and exchange the rows (e.g., Anderson et al., 1991). An example is given in Fig. 5.
Fig. 3. Chromosome decoding when there is no restriction on the delivery of items.
Fig. 4. Chromosome decoding when item 1 can not be shipped with item 2, and item 4 can not be shipped with item 5.
Fig. 5. An example of horizontal crossover.
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
147
Table 1 Values of the parameters.
Fig. 6. An example of horizontal mutation.
(iii) Mutation: In order to prevent the local optimum, the offspring from crossover will be mutated with a certain probability. We apply horizontal mutation where a row of the chromosome will randomly be selected, and then two genes will be swapped randomly (Beasley et al., 1993; Kruse et al., 2013). An example is given in Fig. 6.
Qi ¼
D i mi T ; ð1 E½pi Þ
8 i ¼ 1; …; n:
ð37Þ
Note that when substituting pi ¼ 0; 8 i, the case with no defective items can be analyzed.
4. Numerical examples In order to illustrate the developed model numerically, a JRP is considered with 5 items in which item 1 must not be shipped with item 2, and item 4 must not be shipped with item 5; and there are two trucks each with capacity of 5000 items. The available budget, B, is $40,000 and the ordering cost, K, is $10 per order. The percentage of defective of item i, pi, is assumed to be uniformly distributed, i.e., pi U½vi ; wi , hence E½pi ¼ ðvi þ wi Þ=2 and Var½pi ¼ ðwi vi Þ2 =12. Moreover, the values of the other parameters are given in Table 1: The problem is solved by GA and the results are summarized in Table 2. Items 2 and 5 are assigned to truck 1, items 1 and 4 are assigned to truck 2, 35% of the order quantity of item 3 is assigned to truck 1 and the rest is assigned to truck 2. The constraint that item 1 must not be shipped with item 2, and item 4 must not be shipped with item 5 has been satisfied. Items 1 and 3 are ordered
1
2
3
4
5
Di (units/year) ci ($/unit) μi ($/unit) νi ($/unit) ai ($/delivery) hi ($/unit/year) xi (units/year) E½pi
1000 6 0.15 3 30 1.2 1500 0.01
4000 2 0.1 1 20 0.8 6500 0.02
3000 5 0.15 2.5 80 1 5000 0.015
5000 8 0.2 4 40 3.2 7000 0.04
2500 4 0.15 2 20 0.8 4000 0.030
Table 2 Summary of the results.
3.4. Initialization The initialization process randomly generates the first generation of chromosome equal to the number of population. The generated chromosomes will be decoded so that the solutions and fitness functions are determined. Next, the stopping criterion is checked. If it is not satisfied, then the new generation will be generated by two branches. During this process, first, the GA will preserve some good solutions through elitism by cloning the best chromosomes from the previous generation with a certain probability ðERÞ. Second, the parent chromosomes will be selected by two methods: RRW and TS. Third, the crossover will be performed to the selected parent chromosomes and the mutation will be applied ðMRÞ. This implies that some children after crossover may not be mutated, and finally, the total number of offspring will be equal to the number of population ðPÞ. Hence, at the combined arrow into the “new generation” in Fig. 1, P½MRð1 ERÞ þ ð1 MRÞð 1 ERÞ þ ER ¼ P. After that, the generation number is increased by one. The process will be repeated until the stopping criterion is satisfied as shown in Fig. 1. The stopping criterion is that the difference of fitness function between the current and the last generation is less than 5. The optimal value of family cycle time, T, the reorder frequencies, m0i s, the proportion of each item transported by each truck, bij, are obtained. The order quantity of item i is given by the following equation:
Item
Item
1
2
3
4
5
mi bi1 bi2 Yi1 Yi2
3 0 1 0 1
1 1 0 1 0
3 0.35 0.65 1 1
1 0 1 0 1
2 1 0 1 0
every 3T n , items 2 and 4 are ordered every T n , and item 5 is ordered every 2T n , where T n ¼ 0:079146 year. The expected total cost per year is $88,035.30. Now, in order to understand how defective items affect the ordering policy, we change the percentage of defective items for each item. It is assumed that vi ¼ 0; 8 i ¼ 1; 2; …; 5, and wi takes different values for different items, where 0:02 r w1 r 0:26, 0:05 r w2 r 0:29, 0:035 r w3 r 0:275, 0:03 r w4 r0:27, and 0:022 r w5 r0:262. It should be noticed that as wi is changed in Table 3, both mean and standard deviation of the uniform distribution are changed. From the first row of Table 3, it can be seen that when the items are perfect (i.e., vi ¼ wi ¼ 0), the order quantity for each item and the total expected cost are the lowest. As the percentage of defective items increases, the optimal family cycle length, T, becomes shorter indicating more frequent order. Moreover, order quantity of each item increases while mi for each item remains the same; and also when percentage of defective items increases the expected total cost per unit time increases. This implies that, while satisfying the constant demand for each item, in order to compensate for the defective items, the order quantity should be increased and also the items should be ordered more frequently. This leads to a higher cost when compared to the JRP with perfect items. Next, in order to evaluate the performance of GA, besides the total cost, three criteria are used:
Convergence score ðCSÞ: It determines how well the solution converges to the best value for all generations of GA (De Jong, 1980). CS can be computed by using the following equation: CS ¼
1 Gen
PGen
g¼1
E½C g ð .,.,.,. Þ
ming E½C g ð .,.,.,. Þ
;
ð38Þ
where Gen is the number of generations, E½C g ð .,.,.,. Þ obtained from Eq. (11) is the best total expected cost of generation g. CS C 1 means that the average fitness value is not different from the best value that occurred in this generation. Then it implies that there is no sufficient gradient to reach to the global optimum. Improvement score ðISÞ: This measure determines how well the solution was improved since the first generation (Beasley
148
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
Table 3 Sensitivity analysis for items with imperfect quality. w1
w2
w3
w4
w5
m1
m2
m3
m4
m5
T
Q1
Q2
Q3
Q4
Q5
E½C
0 0.020 0.040 0.060 0.080 0.100 0.120 0.140 0.160 0.180 0.200 0.220 0.240 0.260
0 0.050 0.070 0.090 0.110 0.130 0.150 0.170 0.190 0.210 0.230 0.250 0.270 0.290
0 0.035 0.055 0.075 0.095 0.115 0.135 0.155 0.175 0.195 0.215 0.235 0.255 0.275
0 0.030 0.050 0.070 0.090 0.110 0.130 0.150 0.170 0.190 0.210 0.230 0.250 0.270
0 0.022 0.042 0.062 0.082 0.102 0.122 0.142 0.162 0.182 0.202 0.222 0.242 0.262
3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 1 1 1 1 1 1 1 1 1 1 1 1 1
3 3 3 3 3 3 3 3 3 3 3 3 3 3
1 1 1 1 1 1 1 1 1 1 1 1 1 1
2 2 2 2 2 2 2 2 2 2 2 2 2 2
0.080732 0.079891 0.079334 0.078767 0.078189 0.077600 0.077000 0.076389 0.075767 0.075135 0.074492 0.073838 0.073174 0.072499
242.20 242.09 242.86 243.61 244.34 245.05 245.74 246.42 247.07 247.70 248.31 248.89 249.46 250.00
322.93 327.76 328.85 329.91 330.96 331.98 332.97 333.94 334.88 335.80 336.69 337.55 338.38 339.17
726.59 731.83 734.20 736.52 738.79 741.00 743.16 745.26 747.29 749.27 751.18 753.02 754.80 756.51
403.66 405.54 406.84 408.12 409.37 410.58 411.76 412.91 414.03 415.11 416.16 417.16 418.13 419.07
403.66 403.90 405.18 406.43 407.66 408.85 410.01 411.14 412.23 413.28 414.30 415.29 416.23 417.14
84,265.21 86,214.74 87,495.49 88,803.07 90,138.33 91,502.16 92,895.47 94,319.22 95,774.41 97,262.09 98,783.35 100,339.3 101,931.2 103,560.3
et al., 1993); it can be computed as follows: IS ¼
E½C 1 ð .,.,.,. Þ ; ming E½C g ð .,.,.,. Þ
ð39Þ
where E½C 1 ð .,.,.,. Þ is the minimum total expected cost of generation 1. A larger value of IS is an indication of a better algorithm. Computational time (CT): We measure it by CPU time from the beginning of GA until a termination condition is reached.
As shown in Table 4, we performed the factorial design with four factors, such as population size (P), mutation rate (MR), elitism rate (ER), and selection method, to determine whether they affect the quality of solutions. The details of the factorial design are given in Appendix A. From Table 5, the overall observation is that different population sizes will result in significantly different outputs for the total expected cost, improvement score, convergence score, computational time, and the number of generations. Observations can also be made for an individual response. For a smaller n, the only parameter that affects the total expected cost is the population size. However, for a larger n, there are more factors that affect the total expected cost, and hence these parameters factors should be varied to obtain the optimal solution. For example, for n ¼ 50, the population size, elitism rate, mutation rate, and selection methods should be varied. In addition, it is observed different replications will lead to different computational time requirements. Finally, for a smaller n, the generation numbers of different population sizes are slightly different, but for a larger n, the generation numbers are the same for all factors. We can use four performance indicators (E½C, IS, CS, CT) to compare the GA efficiency. The observations help us analyze the best set of parameters for GA based on the minimum total expected cost in all scenarios, and the best sets of parameters are P ¼30, ER ¼0.4, and MR ¼0.2. From the results of the first set of experiments, it is observed that for a smaller number of items, n, there is no significant difference between the total expected costs obtained from the two selection methods for different mutation rates, elitism rates, and number in population. However, for a larger number of items, n, there is difference in the total expected cost when different mutation rates, elitism rates, and number in population are used, as given in Table 6, where E½CR is the cost determined when RRW is used and E½CT is the cost determined when TS is used. In addition, it is observed that the second method (tournament selection) provides a slightly lower total expected cost. The second set of experiments is implemented for the selected GA parameters given in Table 7, population size, P ¼30, mutation rate, ER ¼0.4, and elitism rate, MR ¼0.2. There are total of
Table 4 Selected values for parameters and selection methods for the first set of experiments. Parameters
Values
Scenarios ðn; KÞ The number of replications Population size (P) Mutation rate (ER) Elitism rate (MR) Selection method Output
(10, 10), (10, 40), (50, 10), (50, 40) 3 10, 20, 30 0.1, 0.2, 0.3, 0.4 0.01, 0.05, 0.1, 0.2, 0.3 RRW, TS E½C, IS, CT, CS, Gen
5 4 2 ¼40 experiments. For both selection methods (RRW and TS), total expected cost, E½C, the improvement score, IS, the CPU time, CT, and the convergence score, CS, are determined. For example, convergence score for RRW and TS are denoted as CSR and CST , respectively, and the results are given in Table 8. With respect to the total expected cost, it is not significantly different in most scenarios. It seems that TS is slightly better than RRW when n increases. As for the improvement score, the larger value the better. It can be seen that, for a smaller value of n, RRW outperforms. However, when n increases, TS outperforms due to a smaller ratio of ISR =IST . In terms of CPU time, when n is small, RRW requires less computational time than that of TS. On the other hand, when n is equal to 50, RRW requires more computational time than that of TS. With respect to the convergence score, the smaller the better. However, observing the ratios of the convergence scores in Table 8, it is inconclusive to determine which method outperforms. In summary, four different measures are used to evaluate the performance of GA. For a smaller number of products, rank roulette wheel performs better based on the total expected cost, improvement score, and CPU time. However, for a larger number of items, tournament selection performs better. Unfortunately, convergence score does not provide any clear conclusion as to which selection method performs better. This indicates that one has to pay attention to the selection method based on the number of items, n, when GA is used to solve the JRP considered in this research. 4.1. Comparison between genetic algorithm and differential evolution algorithm The heuristic approach, differential evolution (DE) algorithm originally presented in Storn and Prince (1997), has several advantages: (a) it can handle different functions such as non-differentiable, nonlinear, and multimodal; (b) it can obtain promising results in a reasonable amount of time by using a vector
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
149
Table 5 Summary of significant p-value of factor to responses of scenarios. Response p-Value of significant factor for each scenario of ðn; KÞ (10, 10)
(10, 40)
(50, 10)
(50, 40)
E½C
Population (0.031)
Population (0.003)
IS CT
Population (0.003) Replication (0.080), Population (0.000), Selection (0.000). Population (0.000) Population (0.047)
Population (0.000) Replication (0.001), Population (0.000), Selection (0.000). Population (0.000) Population (0.002)
Population (0.000), MR (0.052), Selection (0.000). Population (0.000) Replication (0.000), Population (0.000), ER (0.000), Selection (0.000). Population (0.000), Selection (0.000). –
Population (0.000), ER (0.000), Selection (0.000). Population (0.000) Replication (0.000), Population (0.000), MR (0.000), Selection (0.000). Population (0.000), Selection (0.000). –
CS Gen
Table 6 Comparison of the total cost when n¼ 50 and K ¼ $40. ^ER
0.1
0.2
0.3
0.4
0.1
0.2
0.3
0.4
MR
E½CR
0.01 0.05 0.1 0.2 0.3
1,092,778.9 1,093,383.5 1,092,866.0 1,093,100.1 1,093,062.5
1,093,158.4 1,093,229.5 1,092,908.6 1,092,970.1 1,092,834.4
1,093,121.3 1,092,905.4 1,093,300.0 1,092,869.9 1,092,980.8
1,093,051.0 1,092,836.8 1,093,094.8 1,092,896.0 1,092,795.9
1.000044 1.000418 0.999856 1.000307 1.000266
1.000265 1.000513 0.999959 1.000186 0.999889
1.000159 1.000171 1.000410 1.000040 1.000117
1.000154 0.999772 1.000246 1.000185 0.999783
20
0.01 0.05 0.1 0.2 0.3
1,092,734.4 1,092,858.0 1,092,751.6 1,092,702.1 1,092,691.6
1,093,002.3 1,092,733.1 1,092,682.0 1,092,692.3 1,092,757.2
1,092,810.7 1,092,746.3 1,092,900.7 1,092,860.5 1,092,787.2
1,092,807.2 1,092,678.4 1,092,711.9 1,092,698.4 1,092,837.2
0.999932 1.000096 1.000032 0.999950 0.999994
1.000302 1.000067 0.999939 0.999990 1.000047
1.000119 0.999970 1.000199 1.000156 1.000090
1.000088 0.999992 1.000054 1.000015 1.000079
30
0.01 0.05 0.1 0.2 0.3
1,092,898.7 1,092,717.9 1,092,667.9 1,092,665.3 1,092,774.8
1,092,945.3 1,092,742.3 1,092,663.4 1,092,739.8 1,092,700.9
1,093,025.3 1,092,771.8 1,092,659.1 1,092,725.7 1,092,720.9
1,092,864.1 1,092,702.6 1,092,735.1 1,092,684.0 1,092,699.8
1.000200 1.000060 1.000014 1.000012 1.000061
1.000245 1.000070 0.999981 1.000054 1.000030
1.000287 1.000081 1.000004 1.000039 1.000037
1.000176 1.000035 1.000072 0.999972 0.999948
P 10
E½CR =E½CT
population; (c) it has the self-organizing scheme that takes the difference vector of two randomly chosen population vectors to perturb an existing vector; and (d) it has the good convergence properties. The framework of DE is given in Fig. 7. Each step of the framework is described below. 1. Initialization: DE utilizes multi-dimensional parameter vectors to conduct the parallel direct search. Hence, we developed the structure of parameter vectors for two decision variables (mi and bij) as follows: 0
m1
B B m2 B B ⋮ @ mn
b11
b12
⋯
b21
b22
⋯
⋮
⋮
⋯
bn1
bn2
⋯
b1J
Parameters
Values
n K Selection method
5, 10, 20, 30, 50 10, 20, 30, 40 RRW, TS
Eq. (41): ( FðgÞ ¼ 2λ F 0
1
C b2J C C: ⋮ C A bnJ
ð41Þ
λ ¼ e1 Gen þ 1 g ; Gen
Then, we randomly generated the initial values of the vectors according to the uniform distribution U[0, 1] and decoded the values of mi, bij, and Yij as in Section 3.1. 2. Mutation: DE usually adds the scaled difference of two randomly selected parameter vectors to generate a new parameter vector as follows:
νk ðg þ 1Þ ¼ xr1 ðgÞ þF½xr2 ðgÞ xr3 ðgÞ;
Table 7 Selected values for parameters and selection methods for the second set of experiments.
i a r 1 a r 2 ar 3 ;
ð40Þ
where νk ðg þ 1Þ denotes kth new vector for the (g þ1)th generation; and xr1 ðgÞ, xr2 ðgÞ, and xr3 ðgÞ represent r1th, r2th, and r3th target vectors of gth generation, respectively. F is a scale factor that controls the amplification of the differential variation of two target vectors. The decision variables, mi and bij, are mutated based on Eq. (40). The scale factor F is used based on
where F(g) is the scale factor of gth generation, and Gen is the total number of generation. 3. Crossover: In order to increase the diversity of the vectors, crossover generates a trial vector as shown in Eq. (42): uk ðg þ1Þ ¼ ðu1k ðg þ 1Þ; u2k ðg þ 1Þ; …; uzk ðg þ 1ÞÞ; where uzk ðg þ 1Þ ¼
(
νzk ðg þ 1Þ if rand r CR; xzk ðgÞ
if rand 4 CR:
;
z ¼ 1; 2; …; Z;
ð42Þ
ð43Þ
where rand is a uniform random number within the interval [0,1], CR is the user-defined crossover rate, Z denotes the number of dimension of vectors. The decision variables, mi and bij, are crossovered according to Eqs. (42) and (43). 4. Selection: In this step, the trial vector will be compared with the target vector based on the total expected cost function f. The
150
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
Table 8 Comparison of performance measures and for selected GA parameters. n
E½CR
ISR
CTR
CSR
GenR
E½CR =E½CT
ISR =IST
CT R =CT T
CSR =CST
GenR =GenT
5 10 20 30 50
88,035.3 206,909.5 446,537.9 645,911.6 1,092,277.0
1 1.000114 1.000134 1.000875 1.003289
0.5332 0.111575 1.034626 2.181036 8.015088
1 1.000022 1.000011 1.000098 1.000720
10 11 61 81 101
1 1 1 1 1.00004
1 0.999967 0.999934 1.000118 0.997390
4.299147 0.723527 0.669337 0.728164 1.048020
1 1.000006 0.999996 1.000033 0.999893
1 1.1 1 1 1
20
5 10 20 30 50
88,159.1 207,027.0 446,656.8 646,061.5 1,092,402.0
1 1.000079 1.000233 1.000648 1.006614
0.082572 0.106939 1.056045 2.196840 8.045402
1 1.000017 1.000014 1.000092 1.001205
10 10 61 81 101
1 1 1 1.000014 1.000028
1 1.000079 0.999935 1.000015 1.002205
0.668619 0.677760 0.687305 0.735720 1.049723
1 1.000017 0.999999 1.000045 1.000536
1 1 1 1 1
30
5 10 20 30 50
88,278.01 207,142.0 446,774.5 646,182.6 1,092,573.0
1 1.000085 1.00027 1.001048 1.005727
0.082319 0.120387 1.052576 2.194502 8.072410
1 1.000018 1.000015 1.000111 1.000898
10 10 61 81 101
1 1 1 1 1.000038
1 1 1.000091 1.000025 0.999572
0.661697 0.759408 0.684184 0.733724 1.053705
1 1.000004 0.999994 1.000051 0.999889
1 1 1 1 1
40
5 10 20 30 50
88,375.6 207,254.8 446,891.0 646,301.3 1,092,653.0
1 1.000002 1.000152 1.000578 1.006092
0.084821 0.119709 1.034572 2.199863 8.007101
1 1 1.000025 1.000078 1.000812
10 10 61 81 101
1 0.999999 1 0.999956 0.999968
1 0.999953 0.999971 0.999972 1.001295
0.685297 0.748792 0.667291 0.730768 1.056270
1 0.999988 1.000013 1.000025 1.000240
1 1 1 1 1
K 10
be seen that, for a larger value of n, GA obtained a lower total expected cost but required more computational time than that of DE. Nonetheless, overall, the total expected cost obtained from GA is 2.31% less than that of DE. Similarly, computational time (CT) required by GA is 0.52% less than that of DE. Hence, GA is better in terms of these two measures. IS and CS scores determine how well the solutions are improved in subsequent generations. Table 9 shows that most IS values obtained from DE are larger than those obtained from GA, and hence DE is better. Similarly, CS scores obtained from DE are close to one and it implies that the average value of fitness functions is not different from the best value. In addition, most CS scores obtained from DE are slightly larger than those obtained from GA, and hence DE is better. After all, the main objective of the comparison is the total expected cost, and hence the GA is preferred.
5. Conclusion Fig. 7. DE algorithm.
vector with the smaller value of the total expected cost will be put into the next generation ( uk ðg þ 1Þ; f ðuk ðg þ1ÞÞ r f ðxk ðgÞÞ; xk ðg þ1Þ ¼ ð44Þ Otherwise: xk ðgÞ The decision variables, mi and bij, used Eq. (44) to select the child vectors of the next generation. 5. Stopping criterion: We used the same stopping criterion as in the GA. The proposed JRP problem was solved using the DE algorithm and then the results obtained from DE and GA with selection method of rank roulette wheel are compared as shown in Table 9. Experiments are implemented for different values of K and n. The recorded outputs are the total expected cost, improvement score, computational time, and the number of generations. In order to compare the results, differences of outputs from both GA and DE are presented in the last five columns of Table 9. For example, the difference in the total expected cost for K ¼ 10 and n ¼5 is calculated as fð88; 035:3–90; 247:18Þ=88; 035:3g 100 ¼ 2:51%. It can
Resource constraints are always inherent with a JRP; there have been several studies incorporating a few of them. This research simultaneously takes into account not only the constraints often addressed in the literature, such as budget constraint and capacity constraint, but also dependency of items and defective items in each shipment. This adds a higher dimension to the problem and captures the nature of a real JRP problem existing in many retail companies. In order to investigate the ordering policy, the optimal family cycle length, the best set of order frequencies, the order quantity, and the total expected cost per unit time are determined. A two-dimensional genetic algorithm (GA) is proposed to solve the problem. An empirical study is conducted to determine the best GA parameters and selection method. The results also show that a higher percentage of defective items leads to a lower optimal family cycle length, a higher order quantity for each item, and a higher total expected cost per unit time. The JRP is also solved by using the differential evolution (DE) algorithm and results obtained from DE are compared with those obtained from GA. Based on the minimum total expected cost of the JRP, GA is preferred. This research can be extended in several ways. There are a number of ways to add more aspects to the investigated problem.
Table 9 The comparison between GA and DE. n
GA-RRW
DE
Difference (%)
E½C
IS
CT
CS
Gen
E½C
IS
CT
CS
Gen
E½C
IS
CT
CS
Gen
5 10 20 30 50
88,035.3 206,909.5 446,537.9 645,911.6 1,092,277.0
1 1.000114 1.000134 1.000875 1.003289
0.5332 0.111575 1.034626 2.181036 8.015088
1 1.000022 1.000011 1.000098 1.00072
10 11 61 81 101
90,247.18 211,120.9 454,091.3 654,215.5 1,134,355
1 1.01712 1.015416 1.019667 1.039289
0.124276 0.177396 1.553928 2.806978 2.030138
1 1.007874 1.005948 1.005541 1.010408
10 11 61 81 101
2.51 2.04 1.69 1.29 3.85
0.00 1.70 1.53 1.88 3.59
76.69 58.99 50.19 28.70 74.67
0.00 0.79 0.59 0.54 0.97
0 0 0 0 0
20
5 10 20 30 50
88,159.1 207,027.0 446,656.8 646,061.5 1,092,402.0
1 1.000079 1.000233 1.000648 1.006614
0.082572 0.106939 1.056045 2.19684 8.045402
1 1.000017 1.000014 1.000092 1.001205
10 10 61 81 101
90,298.52 211,171.6 453,001.9 654,282.5 1,159,125
1 1.011971 1.011447 1.019677 1.024488
0.121747 0.161616 1.554122 2.828463 2.089074
1 1.005273 1.002913 1.006586 1.005818
10 10 61 81 101
2.43 2.00 1.42 1.27 6.11
0.00 1.19 1.12 1.90 1.78
47.44 51.13 47.16 28.75 74.03
0.00 0.53 0.29 0.65 0.46
0 0 0 0 0
30
5 10 20 30 50
88,278.0 207,142.0 446,774.5 646,182.6 1,092,573.0
1 1.000085 1.00027 1.001048 1.005727
0.082319 0.120387 1.052576 2.194502 8.07241
1 1.000018 1.000015 1.000111 1.000898
10 10 61 81 101
90,349.37 210,157.8 454,202.4 655,691.3 1,139,811
1 1.012068 1.014232 1.032538 1.040375
0.122499 0.163604 1.552599 2.854533 2.031563
1 1.003739 1.002678 1.010285 1.008621
10 10 61 81 101
2.35 1.46 1.66 1.47 4.32
0.00 1.20 1.40 3.15 3.45
48.81 35.90 47.50 30.08 74.83
0.00 0.37 0.27 1.02 0.77
0 0 0 0 0
40
5 10 20 30 50
88,375.7 207,254.8 446,891.0 646,301.3 1,092,653.0
1 1.000002 1.000152 1.000578 1.006092
0.084821 0.119709 1.034572 2.199863 8.007101
1 1 1.000025 1.000078 1.000812
10 10 61 81 101
90,399.72 210,214.1 453,135.5 6,557,58.4 1,133,157
1 1.015391 1.016618 1.035515 1.041107
0.123717 0.159229 1.554821 2.877044 2.031391
1 1.006917 1.006004 1.01086 1.008839 Average
10 10 61 81 101 2.31
2.29 1.43 1.40 1.46 3.71 1.70
0.00 1.54 1.65 3.49 3.48 12.99
45.86 33.01 50.29 30.78 74.63 0.52
0.00 0.69 0.60 1.08 0.80 0
0 0 0 0 0
10
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
K
151
152
P. Ongkunaruk et al. / Int. J. Production Economics 175 (2016) 142–152
For example, incorporating quantity discounts in this particular model could be an extension. Multiple suppliers could also be considered in future research. Moreover, instead of considering complete screening of the shipment, an acceptance sampling approach can be implemented. In terms of solution methodology, differential evolution algorithm can be further improved and then can be compared with the proposed GA. In this research, we considered that the percentage of defective of an item follows a uniform distribution. However, the percentage of defective of an item can follow any other distribution and such consideration can be investigated in the future research.
Acknowledgment The authors would like to thank Natural Sciences and Engineering Research Council of Canada and National Nature Science Foundation of China for funding this research under the Grant number of 71402149.
Appendix A. The details of the factorial design Each factor is varied to the levels that they will more likely affect the results. For example, population size is varied at three levels: 10, 20 and 30. Elitism rate is varied at four levels: 0.1, 0.2, 0.3, and 0.4. Mutation rate is varied at five levels: 0.01, 0.05, 0.1, 0.2, and 0.3, and it is noted that if elitism rate is more than 0.3, the result converges very quickly, and two selection methods are considered: Rank Roulette Wheel (RRW) and Tournament Selection (TS). Then, GA is repeated three times for four combinations of items and the ordering cost, i.e., (n,K) ¼{(10,10), (10,40), (50,10), (50,40)}. The analysis of variance is performed for five responses: the total expected cost ðE½CÞ, improvement score (IS), cpu time (CT), convergence score (CS), and the number of generations (Gen). Hence, for four combinations of ðn; KÞ, there will be 3 3 4 5 2 4 ¼ 1440 scenarios to be solved, but each scenario has five responses and hence the total number of responses is 7200. To be efficient, we performed the main effect analysis and determined which factor affects the response and the results are concluded in Table 5. The main effect analysis is tested with 95% confidence interval (i.e., α ¼0.05) and p-value is given within the parentheses. The p-value indicates if the outputs are significantly different for different levels of the factors.
References Anderson, C., Jones, K.F., Ryan, J., 1991. A two-dimensional genetic algorithm for the ising problem. Complex Syst. 5 (3), 327–333. Beasley, D., Bull, D.R., Martin, R.R., 1993. An overview of genetic algorithms: Part 1, fundamentals. Univ. Comput. 15 (2), 58–69. Bortfelldt, A., 2006. A genetic algorithm for the two-dimensional strip packing problem with rectangular pieces. Eur. J. Oper. Res. 172 (3), 814–837. Bortfelldt, A., Winter, T., 2009. A genetic algorithm for the two-dimensional knapsack problem with rectangular pieces. Int. Trans. Oper. Res. 16 (6), 685–713. Cárdenas-Barrón, L.E., 2000. Observation on: “economic production quantity model for items with imperfect quality” [Int. J. production economics 64 (2000) 59–64]. Int. J. Prod. Econ. 67 (2), 201.
Cha, B.C., Moon, I.K., 2005. The joint replenishment problem with quantity discounts under constant demand. OR Spectr. 27, 569–581. De Jong, K., 1980. Adaptive system design: a genetic approach. IEEE Trans. Syst. Man Cybern. 10 (9), 566–574. Goyal, S.K., 1974. Determination of optimum packaging frequency of items jointly replenished. Manag. Sci. 21 (4), 436–443. Goyal, S.K., 1975. Analysis of joint replenishment inventory systems with resource restriction. Oper. Res. Q. 26 (1), 197–203. Goyal, S.K., Cárdenas-Barrón, L.E., 2002. Note on: economic production quantity model for items with imperfect quality—a practical approach. Int. J. Prod. Econ. 77, 85–87. Hariga, M., 1994. Two new heuristic procedures for the joint replenishment problem. J. Oper. Res. Soc. 45 (4), 463–471. Herbert, E.A., Dowsland, K.A., 1996. A family of genetic algorithms for the pallet loading problem. Ann. Oper. Res. 63 (3), 415–436. Hong, S.P., Kim, Y.H., 2009. A genetic algorithm for joint replenishment based on the exact inventory cost. Comput. Oper. Res. 36, 167–175. Jia, W., Jiang, L., Li, H., Zheng, G., Li, X., 2009. Genetic algorithm for band gap optimization under light line in two-dimensional photonic crystal slab. Opt. Appl. 39 (3), 481–488. Kaspi, M., Rosenblatt, M.J., 1983. An improvement of Silver's algorithm for joint replenishment problem. IIE Trans. 15 (3), 264–267. Kaspi, M., Rosenblatt, M.J., 1991. On the economic ordering quantity for jointly replenished items. Int. J. Prod. Res. 29, 107–114. Khouja, M., Goyal, S.K., 2008. A review of the joint replenishment problem literature: 1989–2005. Eur. J. Oper. Res. 186 (1), 1–16. Khouja, M., Michalewicz, Z., Satoskar, S.S., 2000. A comparison between genetic algorithm and the RAND method for solving the joint replenishment problem. Prod. Plan. Control 11 (6), 556–564. Kruse, R., Borgelt, C., Klawonn, F., Moewes, C., Steinbrecher, M., Held, P., 2013. Computational Intelligence: A Methodological Introduction. Springer-Verlag, London. Maddah, B., Jaber, M.Y., 2008. Economic order quantity for items with imperfect quality: revisited. Int. J. Prod. Econ. 112, 808–815. Moon, I.K., Cha, B.C., 2006. The joint replenishment problem with resource restriction. Eur. J. Oper. Res. 173 (1), 190–198. Olsen, A.L., 2008. Inventory replenishment with interdependent ordering costs: an evolutionary algorithm solution. Int. J. Prod. Econ. 113, 359–369. Ono, T., Ikeda, T., 1998. Optimization of two-dimensional guillotine cutting by genetic algorithms. In: The 6th Proceedings of The European Congress on Intelligent Techniques and Soft Computing, (EUFIT 1998), 7–10 September 1998, Aachen, Germany, pp. 450–454. Papachristos, S., Konstantaras, L., 2006. Economic ordering quantity models for items with imperfect quality. Int. J. Prod. Econ. 100, 148–154. Porteus, E., 1986. Optimal lot sizing, process quality improvement, and setup cost reduction. Oper. Res. 34 (1), 137–144. Robinson, P., Narayanan, A., Sahin, F., 2009. Coordinated deterministic dynamic demand lot-sizing problem: a review of models and algorithms. Omega 37 (1), 3–15. Rosenblatt, M., Lee, H., 1986. Economic production cycles with imperfect production processes. IIE Trans. 18 (1), 48–55. Salameh, M.K., Jaber, M.Y., 2000. Economic production quantity model for items with imperfect quality. Int. J. Prod. Econ. 64, 59–64. Shu, F., 1971. Economic ordering frequency for two items jointly replenished. Manag. Sci. 17, 406–410. Silver, E.A., 1975. Modifying the economic order quantity (EOQ) to handle coordinated replenishment of two or more items. Prod. Inventory Manag. 16, 26–38. Storn, R., Prince, K., 1997. Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces. J. Glob. Optim. 11 (4), 341–359. van Eijis, M.J.G., Heuts, R.M.J., Kleijnen, J.P.C., 1992. Analysis and comparison of two strategies for multi-item inventory systems with joint replenishment costs. Eur. J. Oper. Res. 59 (3), 405–412. Wahab, M.I.M., Jaber, M.Y., 2010. Economic order quantity model for items with imperfect quality, different holding costs, and learning effects: a note. Comput. Ind. Eng. 58 (1), 186–190. Wahab, M.I.M., Mamun, S., Ongkunaruk, P., 2011. EOQ models for a coordinated two-level international supply chain considering imperfect items and environmental impact. Int. J. Prod. Econ. 134, 151–158. Wang, L., He, J., Wu, D., Zeng, Y.-R., 2012. A novel differential evolutional algorithm for joint replenishment problem under interdependence and its application. Int. J. Prod. Econ. 135, 190–198. Wildeman, R.E., Frenk, J.B.G., Dekker, R., 1997. An efficient optimal solution method for the joint replenishment problem. Eur. J. Oper. Res. 99 (2), 433–444.