Optimizing multi-supplier multi-item joint replenishment problem for non-instantaneous deteriorating items with quantity discounts

This paper deals with a new joint replenishment problem, in which a number of non-instantaneous deteriorating items are replenished from several suppliers under different quantity discounts schemes. Involving both joint replenishment decisions and supplier selection decisions makes the problem to be NP-hard. In particular, the consideration of non-instantaneous deterioration makes it more challenging to handle. We first construct a mathematical model integrated with a supplier selection system and a joint replenishment program for non-instantaneous deteriorating items to formulate the problem. Then we develop a novel swarm intelligence optimization algorithm, the Improved Moth-flame Optimization (IMFO) algorithm, to solve the proposed model. The results of several numerical experiments analyses reveal that the IMFO algorithm is an effective algorithm for solving the proposed model in terms of solution quality and searching stableness. Finally, we conduct extensive experiments to further investigate the performance of the proposed model.


Introduction
In real life, the decay or deterioration of products is a natural phenomenon. The majority of sales revenue of supermarkets and grocery stores derives from deteriorating merchandises. For example, deteriorating items contribute 50% of retailing sales in North American and 30% of supermarket sales worldwide. (Chen et al. 2014 [1]). With the development of preservation technologies, retailers could guarantee products not to deteriorate immediately upon produced, resulting in a longer shelf life. This is now commonly applied for almost all food items (e.g., firsthand vegetables, milk, meats and sea foods), pharmaceuticals, fashionable products, and electronic equipment. Wu et al. (2006) [2] initially introduced the term of "non-instantaneous deterioration" to define the phenomenon, and developed an inventory replenishment model for non-instantaneous deteriorating items with stock-dependent demand and partial backlogging. Subsequently, a large number of studies have emerged to consider the inventory problems for non-instantaneous deteriorating items under various settings. However, seldom researchers consider the multi-item inventory policy for non-instantaneous deteriorating The remainder of the paper is organized as follows. Section 2 reviews the related literature. In Section 3, assumptions, notations and the MS-JRNID model formulation are presented. Section 4 introduces the original MFO algorithm and the proposed IMFO algorithm. Numerical examples and parameter sensitive analyses are conducted in Section 5. Section 6 draws the conclusion and depicts the future research.

Literature review
Our study is of relevance to three research streams: inventory models of non-instantaneous deteriorating products, the research of the JRP model as well as its extensions and heuristics to solve these models, and the applications of the MFO algorithm.

Inventory models for non-instantaneous deteriorating items
Given the practical and theoretical relevance, inventory models of non-instantaneous deteriorating products have been extensively investigated since the work of Wu et al. (2006) [2]. As pioneer researchers, they proposed an inventory replenishment model for non-instantaneous deteriorating items with stock-dependent demand to determine the optimal replenishment policy. Thereafter, scholars have focused on inventory problems of the non-instantaneous deterioration items under varies conditions. Tat et al. (2015) [10] proposed an economic order quantity model for non-instantaneous deteriorating items with and without shortages to investigate the performance of a vendor-managed inventory system. Ghasemi (2015) [11] developed a classical economic production quantity model for non-instantaneous deteriorating items by considering a relationship between the holding cost and the ordering cycle length. Pal et al. (2018) [12] developed an inventory model for non-instantaneous deteriorating items with a preservation technology and a random staring time of deterioration. In their model, shortages are allowed or partially backlogged. Mishra et al. (2020) [13] developed a sustainable supply chain inventory policy with controllable non-instantaneous deterioration and environmental policy for greenhouse items.
All of the above studies focus on single-item problems, while only a few researchers pay attentions to multi-item inventory problems with non-instantaneous deterioration. For example, Li et al. (2007) [14] proposed EOQ-based models of perishable products to investigate the impact of the postponement strategy on retailers. Xu and Xiao (2013) [15] developed a JRP model for non-instantaneous deteriorating items with quantity discounts in a supply chain of one-supplier and one-retailer. With full backlogging allowed, Ai et al. (2017) [9] designed an optimal joint replenishment policy for multiple non-instantaneous deteriorating items under constant demand rate. Our research follows the stream by not only considering the inventory replenishment problem for multiple non-instantaneous deteriorating items, but also in a multi-supplier setting.

JRP models and algorithms
Over the last few decades, the JRP has been extensively investigated. Reviews on related works up to the late 1980s were conducted by Goyal and Satir (1989) [16], as well as Aksoy and Erenguc (1988) [17]. Later on, a comprehensive literature review of the JRP related studies from 1989 to 2005 is available in Khouja and Goyal (2008) [3]. Over the past few years, researchers have focused on the studies on the JRP and its extensions, in which the quantity discount effect has been extensively considered. For example, Cha and Moon (2005) [18] studied a joint replenishment problem with quantity discount and developed two efficient heuristic algorithms to solve it. Choudhary and Shankar (2013) [19] developed an integrated model considering all-unit quantity discount for inventory lot-sizing, supplier selection, and carrier selection problem. Paul et al. (2014) [20] developed the JRP models for defective items with or without quantity discount, and proposed an effective heuristic algorithm to solve the models. Ongkunarul et al. (2016) [21] proposed a JRP model with resource constraints and defective items, and developed a genetic algorithm (GA) to determine the optimal reordering policy. Liu et al. (2018) [22] proposed a constrained joint replenishment and delivery model considering quantity discount, and designed a heuristic and hybrid Tabu search algorithm to solve it. Cui (2018) [23] proposed a novel multi-item joint replenishment problem with multiple-type discounts and developed a heuristic algorithm to find the optimal solution of the JRP model. Following these studies, our research focuses on a new JRP model, in which multiple noninstantaneous deteriorating items are joint replenished from multiple suppliers with different quantity discounts. Specifically, although the supplier selection problem is NP-hard, Chakraborty et al. (2020) [24] and Badi and Pamucar (2020) [25] proposed an excellent algorithm to address this complicated problem.
As various JRP models have received considerable critical attentions, many heuristic algorithms have been developed to solve them, such as the RAND algorithm (Cha and Moon, 2005 [18]), the hybrid genetic algorithm (HGA) (Moon et al., 2008 [5]), the differential evolution (DE) algorithm (Wang et al., 2012 [26]), the particle swarm optimization (PSO) algorithm (Kang et al., 2017 [27]), and the modified hybrid Tabu search (TS) algorithm (Liu et al., 2018 [22]), etc. Similarly, this paper follows this stream by developing a new and effective heuristic algorithm to solve the new proposed model.

The MFO algorithm and its applications
Recently, a novel swarm intelligence algorithm called the moth-flame optimization (MFO) has been proposed (Mirjalili, 2015 [7]). It has been then widely applied in different kinds of complicated optimization problems by many researchers, since its high efficiency and low memory usage endows the superiority over other algorithms. For example, Allam et al. (2016) [28] exploited the MFO algorithm in the parameter extraction process of the three diode models for the multi-crystalline solar module. It is noted that the MFO algorithm converges quickly but is easy to fall into a local optimum (Zhang et al., 2016 [8]). Savsani and Tawhid (2017) [29] designed an effective Non-dominated Sorting Moth-flame Optimization (NS-MFO) algorithm that is able to generate non-dominated solutions and true Parato front. They found that the NS-MFO outperforms both the GA and DE algorithm through numerical experiments.  [33] proposed an improved MFO algorithm with adaptive Levy-Flight perturbations to minimize the specific fuel consumption of fighters. Inspired by the effectiveness and vast applications of the MFO algorithm, our study follows this stream by proposing an IMFO algorithm to extend its applications to a new JRP problem.
From the above literature review, the main contributions of this paper are as follows. (a) A new multi-supplier multi-item joint replenishment model for non-instantaneous deteriorating items with quantity discount is proposed. (b) An improved MFO algorithm is designed to solve the proposed model, obtaining results with satisfactory accuracy and robustness. (c) Effects of the key parameters are analyzed to help managers make better decisions on dealing with the selection of several suppliers and the joint replenishment of perishable products.

Assumptions and notations
A MS-JRNID problem involves the coordinated replenishment for several heterogeneous noninstantaneous deteriorating products under different quantity discount. The objective is to find a policy to minimize the total cost. Fig 1 presents a simplified example for illustrating the MS-JRNID network containing three suppliers and four perishable products in a supply chain. In the joint replenishment problem, a major ordering cost will be incurred whenever an order is placed. It is independent of the number of different items in the replenishment order. While a minor ordering cost, involving the packing charge, the preservation cost or order-processing cost, is incurred for each specific item. Obviously, the minor ordering cost may be different when replenishing from different suppliers.
The following assumptions are considered in our problem.
1. The demand rate of each item is constant and deterministic.
2. The replenishment lead time is zero and shortages are not permitted.
3. t di is the fixed preservation period, in which no deterioration of item i occurs. After this period, item i deteriorates at a deterioration rate θ i , where 0<θ i <1. There is no repair or replacement of the deteriorated items during the replenishment cycle.
4. The price of each item is dependent on the magnitude of the replenishment of each item from each supplier. All-units discount scheme is considered in this paper.
5. Each item can be purchased from only one supplier.
6. Under periodic review policy, all items' replenishment is set on a basic cycle time T, and the replenishment T i of Item i is k i times of this basic cycle time, i.e., T i = k i T.
7. Every T time units, there is an order for at least one item.
The following notations in Table 1 are used to develop the model.

PLOS ONE
Optimizing multi-supplier multi-item JRP for non-instantaneous deteriorating items with quantity discounts

Mathematical model
The non-instantaneous deteriorating items begin to decay after the fixed preservation period t di . Hence, there are two possible cases based on the value of t di and T i : Without loss of generality, we focus on a certain non-instantaneous deteriorating item i. During the time interval [0, t di ], the inventory level decreases only due to the demand depletion. Then the item begins to decay and the inventory level drops to zero due to both the demand depletion and the deterioration during the time interval [t di , T].
Hence, the inventory level of item i could be described as follows: ( With the boundary condition I 1i (T i ) = 0, solving the above equations, we get the solution of Eq (1) as follows: The order quantity for each order is then given by: The average holding cost is: the average deterioration cost is: and the total ordering cost is: Suppose that C 1ij , which is a step function of T and k i , is the unit purchasing cost function of item i from jth supplier. For the all-units quantity discount scheme, C 1ij is represented as follows: where y represents the price discount categories. The preceding analysis shows that P ij1 > P ij2 > � � � > P ijy 0 > � � � > P ijy . Furthermore, the marginal growth of quantity discount decreases normally. Therefore, we can obtain the following inequality: As a result, the average purchase cost is denoted by: Therefore, the total cost is: In this case, the product replenishment cycle is less than or equal to the deterioration free time of the product, implying that the retailer can sell out all products before they start to decay. Therefore, we do not need to consider the deterioration. The model turns into the traditional inventory model and the total cost per unit time is given by: where, Combining Case 1 and Case 2, along with T i = k i T, the average total cost is: where k i 2 N n and k i ¼ 1 for some 1 � i � n;

The original Moth-Flame Optimization (MFO) algorithm
The MFO is a novel swarm intelligence optimization algorithm originally presented by Mirjalili (2015) [6]. It is a nature-inspired optimization technology that simulates the transverse orientation mechanism of moths.
In the transvers orientation, the moth flight process is performed with a fixed angle to the moon. Because the moon is far away from the moth, the transverse orientation guarantees flying in a straight line. In reality, moths usually fly spirally around the artificial light. Because once moths see the artificial light, they will try to maintain a similar angle with this light to fly in a straight line. However, such light is very close with respect to the moon, thus moths fly in a spiral path. Moths are deceived by artificial light due to the transverse orientation inefficiency. The transverse orientation inefficiency confirms that it is available only for flying a long distance in a straight path.
The moths have spiral flight around the artificial light keeping similar angle with it. The moths' positions are updated with their spiral flights. And in the optimization process, the spiral flight space must be within the search range. In the MFO, decision variables to be evaluated are the moths' positions in the space. The solutions are the moths whose flights may be one, two, three or hybrid dimensions. Based on a swarm algorithm, the moths' positions can be represented by the matrix as follows: where n is the number of moths, and d is the problem dimensions. Moreover, Eq (11) denotes positions of the moths. Then the corresponding fitness functions are calculated. The vector OM stores the corresponding fitness values of n moths.
Another key component of the MFO is the flames positions defined as: It is shown that the moths and the flames have the same search space dimensions. As for the flames, the fitness values are represented as follows: Moths and flames are both solutions, but are different in the way of being treated and updated. Moths represent the bodies moves in the search space, while flames represent the best positions acquired at present. Furthermore, the moths search and update the variables to be optimized, ensuring that moths will never miss their optimal solution. The MFO can be modelled in three main functions as follows: where: I is a function responsible for both generating random population of moths and calculating the corresponding fitness function; P is the main function of the optimizer that receives the matrix of moths and gives its updated one; and K is a function that verifies the criterion achievement, either true or false. The moth's position is updated according to the flame as follows: where M i is the ith moth, F j is the jth flame and S is the spiral function.
In the original MFO, the moth's position is updated with the logarithmic spiral function. The characteristics of the spiral flight mechanism are as follows: 1. The starting point of spiral flight is the moth.
2. The endpoint of spiral flight is the flame position.
3. The spiral flight space must be within the search range.
For updating the position of each moth and simulating the flight mode of the moth, the MFO are as follows: Thereby, D i can be expressed as follows: For improving the local mining capacity of the moth in the later iteration, it can decrease the number of flames during each iteration in an adaptive linear way. It can be expressed as follows: where l is the current iteration times, N is the maximum number of flames and T denotes the maximum iteration times.

The improved MFO algorithm
In order to enhance the ability of the MFO and to avoid falling into a local optimum, the improvements of the MFO algorithm are executed in the following two aspects.
(1) The original MFO algorithm applies the logarithmic spiral function to simulate the moth's spiral flight path. However, there exist other spiral functions, such as the Archimedean spiral, the Euler spiral, Hyperbolic spiral, the Fibonacci spiral, and the Fermat's spiral etc. After a large number of trials for comparing the optimization performance of different spiral functions, the optimal spiral function could be chosen to simulate the moth's spiral flight path. The test results show that the MFO algorithm with hyperbolic spiral function has the ability to obtain a set of solutions with good convergence and strong distribution in our problem. Therefore, the optimization accuracy and speed could be improved.
(2) Inspired by the Levy-Flight, we propose a MFO algorithm based on it. Because the Levy-Flight has a fast-growing variance, it can be more effective than the brown random motion in a large-scale searching process. In the Levy-Flight strategy, short-distance walking with small step size and occasional large-step walking are alternate. Therefore, on the one hand, some individuals search near the current optimal position to accelerate the local exploitation; on the other hand, the other part of individuals can search in the space far enough from the current optimal position to avoid all individuals falling into the local optimum. Thus, it can escape from the local optimum and improve the global search capability. Fig 2(A), the reciprocal spiral uses the polar equation r = a/ψ, where, r and ψ are the radius and azimuthal angle in a polar coordinate system, respectively. As ψ increases, the spiral winds around the origin and moves closer to it. It winds faster and faster around as it approaches the pole. Fig 2(B) shows definitions of the sector (light blue) and the polar slope angle. Taking the pole as the center of inversion, the hyperbolic spiral inverts to the spiral of Archimedes.

The hyperbolic spiral function. As shown in
Applying the hyperbolic spiral function to simulate the moth's spiral flight path, the Eq (17) can be improved as follows: where S is the hyperbolic spiral function (movement). t is the random number in range [-1, 1], and defines how much the next position of moth should be close to the flame. t = -1 indicates the closest position to the flame while t = 1 indicates the farthest position.

Levy-Flight algorithm.
The Levy-Flight is a probability distribution proposed by French mathematician Paul Pierre Levy. It is a Markov process, where its distribution follows the heavy probability distribution. Inspired by the Levy-Flight, many scholars applied the Levy-Flight strategy to improve their algorithms, which effectively improve the quality of the solution.
The Levy-Flight mechanism can be defined as follows: where s is the random moving step size, and λ is the exponential parameter. Eq (21) is the heavy tail probability distribution, which is hard to performance through the simple programming language. Therefore, when calculating the Levy-Flight searching distance, it is common to apply the following Eq (22) to simulate the Levy-Flight path proposed by Mantegna.
where s is the Levy-Flight distance, and parameter β is defined in range (0, 2), and normally set as β = 1.5. Both μ and ν represent normal distribution showing in the Eq (23), and the standard deviation corresponding to the Eq (23) satisfies to the Eq (24).
As a result, the Levy-Flight path can be calculated through Eq (22) to Eq (24). Applying the Eq (25) instead of the Eq (18) in the original MFO algorithm can express the moth's Levy-Flight path as follows: According to the population aggregation status, the Levy-Flight disturbance strategy can adaptively adjust the perturbation probability to enhance the ability of the IMFO algorithm to escape from local optimal solutions and realize the minimum total cost.

The algorithm design
Based on the above discussions, the procedure of our proposed IMFO algorithm for solving the MS-JRNID model is as follows: Step 1: Determine the problem dimensions, the searching scope and the criterions of parameters K i , T, X ij to be optimized.
Step 2: Determine the fitness function of the MFO algorithm. This study chooses Eq (10) as the fitness function.
Step 3: Initialize the parameters. Set the number of iterations, the population size, the search space and the flame number. Set the current iteration times l = 0.
Step 4: Calculate the fitness value of individuals. Find the current optimal moth position and keep it as the flame fitness value matrix. Judge whether the maximum iteration number is reached, if reached, go to Step 8, otherwise, execute Step 5.
Step 5: Execute iterative processes. Update the flame number based on Eq (19), calculate the distance between the moth and the flame, and update the moth position and the flame position based on the Eq (20) and Eq (25).
Step 6: Calculate the fitness value of individuals. Keep the moth position and the flame position based on the Eq (12) and Eq (14).
Step 7: Find the current optimal moth position. If the current position is better than the previous position, keep the current flame position as the optimal position. Judge whether the maximum iteration number is reached or not, if reached, go to the Step 8, otherwise, set l = l + 1 and return to the Step 5.
Step 8: Output the optimal individual and global extremum, that is, the final flame position K i , T, X ij and corresponding fitness value. The process ends.

PLOS ONE
Optimizing multi-supplier multi-item JRP for non-instantaneous deteriorating items with quantity discounts

Numerical experiments
This section aims to evaluate the performance of the proposed approach and MS-JRNID. In order to verify the competitive performance of the IMFO algorithm, it is compared with other five meta-heuristic algorithms, named Genetic Algorithm (GA) (Deb et al. 2002) [34], Grey Wolf Optimizer Algorithm (GWO) (Mirjalili et al. 2014) [35], Fruit-fly Optimization Algorithm (FOA) (Pan, 2012) [36], Particle Swarm Optimization (PSO) Algorithm (Shi, 2001) [37] and original Math-Flame Optimization Algorithm (MFO) (Mirjalili. 2015) [7]. Furthermore, sensitive analyses of the parameters are performed to find the optimal results of the proposed model.
The parameters of the numerical example come from the related works in literature (Moon et al. 2008 [5], Cui et al. 2016 [38], Liu et al. 2018 [22]). Since they did not consider the deterioration of items, we introduce the parameters θ i , t di , and c i , which represent the deterioration rate, the non-deterioration time and the deterioration cost, respectively. The values of parameters are shown in the Table 2. Moreover, the quantity discount strategies provided by each supplier are shown in the Table 3.

Comparison between GA, FOA, PSO, GWO, MFO and IMFO
In this subsection, we compare the performances of the proposed IMFO algorithm with other algorithms. Among the meta-heuristic methods, the GA is a stochastic search algorithm following the mechanism of natural selection and "survival of the fittest", and has been widely used in production and operations management problems (May et al. 2015) [39]. The FOA has been shown to be an efficient and effective algorithm for many JRPs (Wang, Shi, and Liu 2015) [40]. The PSO is a well-known method that optimizes a problem by iteratively trying to improve a candidate solution with regard to a given measure of quality, the experience of which is very convenient for testing other algorithms. The GWO algorithm mimics the leadership hierarchy and hunting mechanism of grey wolves in nature, it is a new meta-heuristic inspired by grey wolves. Therefore, GA, FOA, PSO, GWO and MFO are selected for comparisons with the proposed solution method in this paper.
In order to facilitate a direct comparison of the performance of different algorithms, the control parameters of GA, FOA, PSO, GWO, MFO are selected based on pilot experiments from Deb et al. (2002) [34]; Pan (2012) [36]; Shi (2001) [37]; Mirjalili et al. (2014) [35]; and Mirjalili (2015) [7] respectively. For all algorithms, the population size is kept to be 30. The iteration termination condition is to stop after 500 generations. Each algorithm is run for 20 times with different seeds. All experiments are performed by MATLAB in the computer with CPU@2.27GHz and 2.00GB RAM.

PLOS ONE
Optimizing multi-supplier multi-item JRP for non-instantaneous deteriorating items with quantity discounts Under the given values of parameters, we can obtain the computational results of different algorithms shown in Table 4. The proposed IMFO can find a better solution compared to other algorithms. To display the convergence features of each algorithm in searching processes, the average best-so-far solution progress over iterations number is shown in Fig 4. Based on the computation results shown in Table 4 and Fig 4, the following results can be derived.
1. The basic cycle time T and the optimal replenishment strategy k i show different features under different algorithms.
2. The optimal total cost of the IMFO algorithm is smaller than that of GA, FOA, GWO, PSO, and MFO. The average percentage savings in total costs provided by the IMFO are 0.01-1.35 when compared to other five algorithms.
3. As for the convergence features of each algorithm, it can be seen that GA, PSO and GWO converge quickly in the early stage of the search process, but converge slowly in the late stage of the search process, while the IMFO finds the minimum total cost at the end of iterations. Therefore, the IMFO shows a competitive advantage over other algorithms.

Sensitivity analyses
In this subsection, we explore some managerial insights based on sensitivity analyses of the parameters. We investigate the effects of changes in system parameters D i and h i on the 1. When the values of parameter D i , h i and S increases, the optimal total costs of each algorithm will also increase. In addition, when D i changes, the minimum total cost of IMFO are slightly less than the cost obtained from GWO, PSO, and MFO, and much less than those from GA and FOA.
2. With the changes of h i , the minimum and average total cost obtained from FOA is the highest, while the solution of the IMFO is the best.
3. It can be found that when S changes, as for the results performance of the average total cost, the IMFO is more competitive.
Furthermore, when D i , h i , and S change, the optimal supplier mixture strategy (X ij ), the optimal replenishment strategy (k i � ), the optimal basic replenishment cycle (T � ), and the optimal system total cost (TC � ) are shown in Tables 5-7.
Based on the computational results shown in Tables 5-7, the following observations can be derived: 1. The optimal total cost is weakly positively sensitive to changes of h i and S, whereas the optimal total cost is highly positively sensitive to changes of D i . The demand rate has the strongest effect on the optimal total cost.
2. It can be found that the optimal solution (k i � , T � , X ij � ) does not change regularly with the system parameters. The main reason lies in the fact that the proposed model considers the multi-supplier's mixture strategy, where the different quantity discounts and the different minor ordering cost are obtained from different suppliers.
3. The optimal replenishment strategy (k i � ) and the basic replenishment cycle (T � ) are highly positively sensitive to the change of system parameters, while the optimal supplier-selection strategy remains relatively robust to the change. This indicates that the retailers may only need to redesign T � and k i � , and to keep the optimal solution of X ij � unchanged. This would save a great number of computational efforts while still achieving the optimal performance.

Conclusions
In this paper, we provide a pioneering focus on the JRP for multiple non-instantaneous deteriorating items with quantity discounts offered by multiple suppliers. A multi-supplier multiitem joint replenishment model is developed, aiming to minimize the average total cost. As this problem is NP-hard, a novel meta-heuristic algorithm, i.e., the MFO, is introduced and redesigned to solve the proposed model. From the results of numerical experiments, we can conclude that our proposed algorithm outperforms most well-known algorithms including GA, FOA, GWO, PSO, and MFO in both the solving efficiency and accuracy. Therefore, it could be considered as an alternative to solve optimization problems among the current algorithms. Sensitivity analyses of the key parameters indicate that the optimal replenishment policy and the basic replenishment cycle are highly positively sensitive to the changes of system parameters, while the supplier selection strategy is strongly robust to them.
This study provides a practical approach to daily inventory problems faced by managers who have to frequently replenish deteriorating items from multiple suppliers with quantity discounts. In real life, the coordinated replenishment with one-supplier offering different quantity discounts for different items could be regarded as a special case of the proposed model. To the best of our knowledge, our work is the first to investigate the joint replenishment problem for non-instantaneous deteriorating items with several suppliers offering different quantity discounts. Nevertheless, there are still some limitations in this study. For example, the demand and deterioration functions are assumed to be constant, and these assumptions may restrict the applications of the model. Moreover, the flame individual only records the historical optimal solution of its corresponding moth individual, and the update of the moth individual is always based on its corresponding flame individual. In other words, there is no communication among the moth individuals.
Our study provides several directions for future scholars in this research stream. For example, it would be interesting to investigate the optimal policy with more general demand rates. In addition, the promotional efforts and preservation technology investment for perishable products could be further considered. Furthermore, the algorithm can be extended to solve the JRP with resource constraints (such as carbon emission constraint, limited capital budget, etc.). In future research, once the better algorithms are developed, our IMFO algorithm could act as a valuable benchmarking algorithm.   Table 5. When Di changes, the optimal solution obtained by the IMFO algorithm.
ΔD i (%) The optimal supplier mixture strategy The optimal replenish strategy (k i � ) The optimal basic replenish cycle (T � ) The optimal total cost (TC � )

PLOS ONE
Optimizing multi-supplier multi-item JRP for non-instantaneous deteriorating items with quantity discounts Supporting information S1 Data. The data set. (DOCX)

Author Contributions
Conceptualization: Xueyi Ai. Table 7. When S changes, the optimal solution obtained by the IMFO algorithm.

S
The optimal supplier mixture strategy The optimal replenish strategy (k i � ) The optimal basic replenish cycle (T � ) The optimal total cost (TC � )  Table 6. When D i changes, the optimal solution obtained by the IMFO algorithm.
Δh i (%) The optimal supplier mixture strategy The optimal replenish strategy (k i � ) The optimal basic replenish cycle (T � ) The optimal total cost (TC � )

PLOS ONE
Optimizing multi-supplier multi-item JRP for non-instantaneous deteriorating items with quantity discounts