In the spare parts supply network, there are many uncertain factories, such as unpredictable demands and changeable lead times. The spare parts shortage caused by those uncertainties may lead to severe losses. To solve the uncertainty of supply network, a determined optimization model is developed and then reformulated as a robust counterpart. In the robust model, it is only necessary to know the moment information of the uncertain parameters rather than the true probability distribution. The solution obtained by the robust model can satisfy the constraints in the worstcase, that is, feasible for any probability distribution within the moment based ambiguity set. Two moment based robust models are studied in this work. The result of the experiment indicates that the robustness of the robust model is stronger than that of the determined model and chance constraint model, and the effect of safety tolerance on the robustness is revealed by sensitivity analysis. Finally, the second order moment model is verified to be superior to the first order moment model in spare parts supply network optimization.
Keywords: Spare parts supply network, uncertainty, moment based ambiguity set, chance constraint optimization, distributed robust optimization
Spare parts supply network optimization is one of the hot issues of spare parts management. Many research works have been done on supply network optimization in certain circumstances. In recent years, many researchers took attention to uncertainty optimization. In the spare parts supply network, many factories such as demands, lead times, and risks are uncertain since the difficulty of data collection or complex dynamic scenarios. Although there are many advanced methods, such as reliability based [1] and data based methods [2], have been used to predict the demands or lead times, most of them are based on overideal assumptions and not accurate enough. In fact, in most cases, it almost cannot forecast the demand precisely at all.
To simplify the problems, many researchers assume that the demands are determined [3] or obeys specified known distribution, such as Normal distribution [4], Poisson distribution [5], Stuttering Poisson distribution [6], etc. There are also many works assumed the lead times are fixed [7] or obey known distributions [8]. These assumptions may be not effective sometimes and may cause out of stock or more serious consequences. To handle the uncertainty in the spare parts supply network, most conmen used methods are stochastic program [9], fuzzy optimization [10], and robust optimization [11]. However, in stochastic optimization, it is usually necessary to know the distribution of uncertain parameters in advance or very timeconsuming by using Monte Carlo simulation. The fuzzy optimization cannot guarantee the feasibility in all time. Thus, the robust optimization was adopted in this paper because it doesn’t rely on the distribution of uncertain parameters and can guarantee the feasibility of the solution in the worst case.
This paper optimizes the supply network without thus assumption to avoid the serious losses caused by small probability events. The demands and lead times are uncertain with unknown distribution, the only things can we get from the historical data is certain known moments or known structural properties. A moment based ambiguity set was established to descript all the possible distribution, which must contain the real distribution of uncertain parameters. The robust model only needs to ensure the worst case in the ambiguity set is satisfied. The main contributions of this paper are as follows. First of all, the uncertain parameters can obey any distribution with known moments in this paper. Secondly, uncertain demand and lead time are considered at the same time. Thirdly, the robust counterparts are formulated based on inequality theory, which is easily understood than duality theory.
The remaining sections of this article are organized as follows. In Section 2, the determined model of spare parts supply network optimization is developed and then reformulated as its moment based robust counterpart. In section 3, the numerical experiment and sensitivity analysis are carried out. In section 4, the results of the robust model are compared with the determined model and chance constraint model, respectively. Finally, we summarize the paper in section 6.
The classical three echelon spare parts supply network consists of supply centers, distribution centers and customers. The ordered spare parts are supported by supply centers and then distributed through the distribution centers to meet the spare parts demand of the customers. In the process of supplying, in order to save the costs, it is often unnecessary to open all of the alternative distribution centers. The research aims to decide which distribution center should be opening and also the number of spare parts flow among the nodes of the spare parts supply network.
Different from the traditional supply network, there are many uncertainty parameters in our spare parts supply network. Lead times and demands are the two most conmen considered uncertainty factories in the supply network. In terms of spare parts demands, the predicted values usually do not agree with the real data especially in some complex circumstances. Although, there are lots of stateoftheart forecast methods are studied and used by researchers, there are still big gaps between the mathematic models and the realworld practice. Similarly, it’s also difficult to ensure the lead times of supply remain unchanged every time, even there may be some accidents that occur in transit. One of the general strategies to handle these uncertainties is assuming these uncertain parameters obey specified distribution. This assumption will not accurate because of the limitation of the historical data. So, in our supply network optimization model, the probability distribution of lead times and demands are even typically unknown. There are only some useful information can be obtained from the samples, such as certain known moments or known structural properties.
The proposed model focus on the problems that optimize the spare parts network with uncertain lead times and demands to meet the customers’ demands within the given hard time window and with the minimum total cost. The uncertainty makes the problem more complicated to be handled. Firstly, the determined model is developed, and then the corresponding robust counterpart problem will be reformulated. The objective function of the model is to minimize the total supply costs, and the constraints of the model include demand fill rate, lead time, captain limitation, flow equilibrium, and so on. The parameters and models are described as follows.
: the amount of supply centers, ;
: the amount of distribution centers, ;
: the amount of customers, ;
: the lead time between supply center and distribution center , which is uncertain;
: the lead time between distribution center and customer , which is uncertain;
: the required maximum lead time.
: unit spare parts transportation costs between supply center and distribution center ;
: unit spare parts transportation costs between distribution center and customer ;
: fixed opening costs of distribution center ;
: unit spare parts inventory cost of distribution center ;
: unit spare parts inventory cost of customer ;
: unit shortage loss of customer ;
: spare parts demand of customer, which is uncertain ;
: the maximum spare parts capacity of distribution center .
Decision variables:
: the amount of spare parts supply from supply center to distribution center ;
: the amount of spare parts supply from distribution center to customer ;
: binary variable, used to label the opening of distribution center . If distribution center is opening, , otherwise, .
Firstly, the determined model is developed. In this model, all parameters seem as determined parameters. The objective function is minimizing the total supply costs, which consists of opening cost, transportation costs, inventory costs, and shortage loss. The objective function is formulated as follows:

(1) 
where represents the total opening costs, the total transportation costs, the total inventory costs, and the total shortage loss. Each costs are calculated as follows:

(2) 

(3) 

(4) 

(5) 
where represents .
The following constraints should be satisfied:
s.t.

(6) 

(7) 

(8) 

(9) 

(10) 

(11) 
where, constraint (6) is a hard time window constraint, which specifies that the spare parts must arrive at the customers before the deadline. The constraints (7) and (8) state that the closed distribution centers should not participate in the supply of spare parts. They also specify that the input of spare parts should not exceed the maximum capacity of distribution centers. Constraint (9) specifies that the demands of spare parts of all the customers must be satisfied. Constraint (10) ensures that the output of spare parts should not exceed the spare parts flow into each distribution center. Constraint (11) states the character of decision variables.
In the determined supply optimization model, there will be two situations of spare parts shortage. One of them is the spare parts cannot arrive in time, and another is the amount of supplied spare parts cannot meet the demand of customers. Both the two kinds of shortage can lead to serious consequences. In the determined model, it is only necessary to find the feasible solutions that satisfy all the constraints of the model. However, in the uncertain model, it is hard to find feasible solutions because of the uncertainty of lead times and demands. We need to reformulate the mathematical model in the uncertain environment and find robust solutions to ensure that the demands of spare parts are satisfied in any case. This section reformulates the determined model as an uncertain model and derives a tractable and efficient deterministic robust counterpart for the uncertain model.
For most supply practices, the supply environment is complex and changeable, and the sample data is limited. It is usually difficult to get the true probability distribution of these parameters. What only can we do is that extract useful information from these limited samples and use the statistical features to modeling and analysis. Despite the unknown distribution of parameters, the moment, such as, mean, variance, covariances, etc., can be obtained easily. The known moment of the sample contains a lot of valuable information about the distribution, so we use a momentbased ambiguity sets describe the uncertain parameters. This ambiguity set consists of all possible probability distributions satisfied the moment condition, and deservedly contains the true distribution of our parameters. If we can find the robust solutions that satisfy the constraints in the worstcase, then the constraints with real distribution can also be met by the robust solutions.
This section focuses on the robust counterpart reformulation in the case of known first order moments and the case of known both the first order and the second order moments. The ambiguity set of these two cases are expressed as follows, respectively:

(12) 

(13) 
where, Eqs. (12) is the ambiguity set with known first order moment, and Eqs. (13) is the ambiguity set with the known first and second order moment. The parameter and represent the mean and variance obtained from historical data, denotes the support of the family of all probability distributions with known moment information and is the probability density function (PDF) of the random variable .
If the function objective consists uncertain parameters, it can be use the expectation of objective to reformulate, or transfer the objective function to constraints [12]. In our spare parts optimization model, the uncertain parameters just exist in constraints. In Eqs. (6) and (9), the lead time and and the spare parts demands are uncertainty with some unknown distribution. It can be formulate as chance constraints to ensure that the constraint must be satisfied under a probability which greater than the specified threshold.
The single chance constraint of Eqs. (6) is as follows:

(14) 
The single chance constraint of Eqs. (9) is as follows:

(15) 
where, represents a specified safety tolerance or violation probability, and all constraints take the same .
In terms of random variable , if the probability density function is known, it is easy to solve the inequality by [13]. However, in our spare parts network, the distribution of lead times and demands are unknown. Even the PDF of the lead times is known, it still be difficult and timeconsuming to calculate the PDF of . Therefore, further derivation of Eqs. (14) and Eqs. (15) is needed. We use the robust chance constraints to formulate inequality (14) and (15) as follows:

(16) 

(17) 
In the case of only the first order moments are known, for an arbitrary distributional random variable . According to the Markov's inequality, there is . Therefore, we can get the robust counterparts of Eqs. (16) and (17) as follows:

(18) 

(19) 
Proof. For inequality (16), there is:

For inequality (17), there is:

In the case of both the first and second order moments are known, for an arbitrary distributional random variable . According to the Cantelli's inequality, there is, . Therefore, we can get the robust counterparts of constraint (16) and (17) as follows:


(21) 
Proof. For inequality (16), there is:

For inequality (17), there is:

In summary, the robust counterpart is derivated by replace the inequality (6) and (9) of the determined model by the inequality (18) and (19) or inequality (20) and (21).
In this case, the three echelon spare parts supply network consists of 2 supply centers, 5 distribution centers, and 4 customers. The spare parts are supported to meet the demands of customers. The first and second moments of the lead times and demands are counted from 50 samples. The relative data are shown in Tables 14. The required maximum lead time is 300 hours, and the safety tolerance is set from 0.1 to 0.9, respectively. The problems were solved by CPLEX and differential evolutionary algorithm.
The optimal solutions of the robust optimization model with known first order moment (referred to as the first order moment model) and with known first and second moment (referred to as the second order moment model) are obtained, their corresponding objective functions and constraints are also calculated and shown in Table 5. In order to analyze the effect of safety tolerance on the robustness of the model, we take multiple values of in the first order moment model and the second order moment model, respectively. Constraint 120 in Table 5 is the value of the left side of the corresponding inequality constraints, that is, if the value is big than zero, the corresponding constraints are violated, otherwise, the constraint is satisfied. The following conclusion can be summed from the table:
Distribution centers 1  Distribution centers 2  Distribution centers 3  Distribution centers 4  Distribution centers 5  

Opening costs  1000  1500  2500  1800  2000 
Inventory costs  12  8<  5  10  10 
Inventory capacity  50  55  80  60  65 
customer1  customer2  customer3  customer4  

Shortage loss  100  120  80  150 
Inventory costs  10  15  12  
Expectation of demands  68  61  57  88 
Variance of demands  9<  11<  7</  8 
Distribution centers 1  Distribution centers 2  Distribution centers 3  Distribution centers 4  Distribution centers 5  

Supply center1  105  104  102  104  110 
Supply center2  94  108  105  98  96 
Customer 1  46  42  59  54  46 
Customer 2  41  67  61  67  41 
Customer 3  49  53  51  49  49 
Customer 4  55  49  60  51  55 
Distribution centers 1  Distribution centers 2  Distribution centers 3  Distribution centers 4  Distribution centers 5  

Supply center1  12/1.6  11/2.5  9/1.8  11/1.1  10/2.3 
Supply center2  7/0.7  15/1.3  15/1.6  7.5/1.8  8/1.1 
Customer 1  5/0.2  5/0.3  7/0.5  3.7/1.2  4/0.8 
Customer 2  4/0.8  5/1.1  3/0.2  4.5/0.6  4/0.3 
Customer 3  3.5/0.3  2.5/0.3  2.9/0.7  1.5/0.4  5/0.8 
Customer 4  2.6/0.2  4/0.4  1.6/0.2  2/0.5  3.5/0.7 
Model  First order moment model  Second order moment model  

Safety tolerance  0.1  0.3  0.5  0.7  0.9  0.1  0.3  0.5  0.7  0.9 
Costs  42029  41679  41168  40958  39480  42634  41172  41120  39459  39392 
Constraint 1  130.3  86.3  23.3  30.2  199.8  1421.134  5796.052  10530.2  10105.74  13000.61 
2  0  0  0  0  0  0  3  2  6  5 
3  0  1  0  0  4  0  1  0  2  0 
4  0  0  0  0  0  4  1  5  0  1 
5  0  0  0  0  0  1  0  0  2  2 
6  1  11  0  0  0  3  4  2  6  6 
7  2  2  1  0  0  0  5  9  7  6 
8  1  5  4  1  4  6  4  0  2  0 
9  0  1  2  1  0  4  2  5  0  2 
10  2  0  0  0  0  8  0  1  3  2 
11  7  11  1  0  0  7  5  6  6  6 
12  59.2  41.6  26.5  15.5  0.4  17.5  126  195.5  82  56.7 
13  55.4  44.5  31  16.2  0.2  2.2  3.1  7  3  2.5 
14  51.8  44.4  26.5  19.2  0.6  0.1  5.9  9  0.7  2.9 
15  77.8  56.2  39  7.5  0.2  12.4  30.7  8.5  98.4  128.8 
16  2  2  1  0  0  0  2  7  1  1 
17  1  4  4  1  0  6  3  0  0  0 
18  0  1  2  1  0  0  1  0  0  1 
19  2  0  0  0  0  7  0  1  1  0 
20  6  0  1  0  0  4  1  4  0  0 
To verify the robustness of the uncertain model, we compare the results of the robust model with the determined model and the opportunity constraint model, respectively. Robust model refers to the first or second order moment model. The determined model refers to the model with the expected value of the uncertain parameter, which is adopted by most decision makers. The chance constraint model refers to the model with the assumption that the distribution of uncertain parameters is known. Because the distribution is unknown, the Monte Carlo method is used to solve the chance constraint model.
In the determined model developed in section 2.3 of this paper, the expectation of lead times时and spare parts demands are adopted as the determined parameters. The comparison with the determined model is taken from two aspects. In the first aspect, the solutions of the first order moment model and the second order moment model (referred to as the robust solution) are brought into the deterministic model to test the feasibility of the robust solution in the determined model.
Brought the robust solutions under different safety tolerances into the determined models, and the values of all constraints are calculated respectively. The constraint value indicates constraint violation, and if the constraint value is less than or equal to 0, it means that the constraint is satisfied. Figure 110 is the histograms of the constraints values in these ten different cases, respectively (two different order moment models with five safety tolerance). Figure 1, taken as an example, shows the constraints values that bring the robust solutions of the first order moment model with a safety tolerance of 0.1 into the determined model. The horizontal axis represents the index of the constraint, and there are a total of 20 constraints in our case. The vertical axis is a constraint value, the yellow bars (labeled as robust model) are the constraint values of the robust model, and the blue bars (labeled as certain model) are the constraint values of determined the model.
In Figure 1, the constraints 1,12,13,14 and 15 of the robust model are not satisfied. In the corresponding determination model, there are only two violated constraints, that is, constraint 13 and constraint 14. It can be seen that although the robust solution is still unfeasible in the deterministic model, it greatly reduces the constraints violate degree. It can be seen that the first order moment model has good robustness, and the same conclusion can be drawn from Figures 2 to Figure 10.
Comparing the results in Figures 1 to Figure 5 (or Figures 6 to Figure 10), it can also be found that the constraint violations of both the determined model and the robust model decrease as the safety tolerance increases, which is consistent with conclusion (2) in section 5.2. Comparing the results of first order moment model with those in second order moment model, we can find that the robust solution of the secondorder moment model makes the constraint violations in determined model much smaller than the robust solution of the firstorder moment model. This conclusion is consistent with conclusion (3) in section 5.2.
Figure 1. Constraint violation comparison between first order moment model and determined model as safety tolerance is 0.1 
Figure 2. Constraint violation comparison between first order moment model and determined model as safety tolerance is 0.3 
Figure 3. Constraint violation comparison between first order moment model and determined model as safety tolerance is 0.5 
Figure 4. Constraint violation comparison between first order moment model and determined model as safety tolerance is 0.7 
Figure 5. Constraint violation comparison between first order moment model and determined model as safety tolerance is 0.9 
Figure 6. Constraint violation comparison between second order moment model and determined model as safety tolerance is 0.1 
Figure 7. Constraint violation comparison between second order moment model and determined model as safety tolerance is 0.3 
Figure 8. Constraint violation comparison between second order moment model and determined model as safety tolerance is 0.5 
Figure 9. Constraint violation comparison between second order moment model and determined model as safety tolerance is 0.7 
Figure 10. Constraint violation comparison between second order moment model and determined model as safety tolerance is 0.9 
In the second aspect, take the feasible solutions obtained from the determined model into robust models to test the feasibility of the determined solutions in the robust model, and the results are shown in Figure 11. It can be seen from the figure that all the determined solutions are infeasible in the robust model. The results indicate that the solutions of the determined model can only guarantee the feasibility in a specific case, but not in the models with uncertain distribution parameters. Generally, the constraints violation degree in the second order moment model is significantly smaller than that in the first order moment model, which indicates that the first order moment is stricter than the second order moment model. With the increase of safety tolerance, the constraint violation degree decreases gradually, which indicates that the smaller the safety tolerance, the more strict of the robustness of the model is.
Figure 11. Constraint violation comparison between robust model and chance constraint model 
In this section, the feasibility of robust solutions in the chance constraint model is verified. Take the robust solutions of first order moment model and second order moment model into the chance constraint model to test the feasibility. Since the uncertain parameters only exist in the constraints (6) and (9) of the spare parts supply model, only these two inequality constraints need to be tested. The corresponding opportunity constraint models are:

(22) 

(23) 
Since unknown the probability distribution of lead times and demands, it cannot be solved by probability density function. Even if the distribution of uncertain parameters is known, the solution of the joint probability density function in Equation (23) is very timeconsuming. So the Monte Carlo method is used to verify the feasibility of the solution
Because the first order moment robust optimization model is not feasible, only the second order moment optimization model is compared and analyzed.
50 samples of uncertain parameters are selected, take the robust solution of the secondorder moment model into these samples, and calculate the times and that inequalities (22) and (23) are satisfied

(24) 

(25) 
where, , if is satisfied. else if, . Then, and are used to represent the value of and . If or were greater than or equal to , the constraints (22) and (23) are satisfied.
The constraint satisfactions of the chance constraint model are shown in Table 6. It can be seen that the constraints (22) and (23) are satisfied in all scenarios. As the safety tolerance decrease, the value of and are larger. This is because when the safety tolerance is small, the robust model is more strict, and the robustness of the robust solution is stronger, which lead to the greater the probability that the constraint is satisfied in chance constraint model.
0.1  0.3  0.5  0.7  0.9  

0.9692  0.8997  0.7938  0.7288  0.6992  
,  0.9269  0.8194  0.8187  0.7938  0.7994 
,  0.9524  0.8219  0.8170  0.7922  0.7223 
,  0.9461  0.9403  0.9396  0.8959  0.8843 
,  0.9289  0.8351  0.7378  0.7298  0.6537 
In this paper, the spare parts network optimization with uncertainty lead times and demands is studied. The determined supply model and its robust counterpart are developed, respectively. In the uncertainty model, the probability of the uncertainties are unknown, we focus on the reformulation of this model with the known first and second order moment. The robust solutions of uncertain models are obtained and the sensitivity analysis with different safety tolerance is curry out. Then the robustness is compared between the robust model and determined model and chance constraint model. The main contributions of this paper are as follows. Firstly, the ambiguity set is used to handle the uncertain probability distributions, and the robust counterpart of the supply model is reformulated to make the model tractable and computable. Secondly, the robust solutions that ensure the feasibility in the worstcase are obtained. These solutions ensure that the uncertain constraints are satisfied under any probability distribution within the moment based ambiguity set. The results of the case study show that the robustness of these solutions are stronger than those in the determined model and chance constraint model. Finally, we verified that the second order moment model is superior to the first moment order model because the second order moment model can make a balance in robustness and flexibility.
[1] Eaves A.H.C.，Kingsman B.G. Forecasting for the ordering and stockholding of spare parts. Journal of the Operational Research Society, 55(4):431437, 2004.
[2] Mao H.L., Gao J.W., Chen X.J., Gao J.D. Demand prediction of the rarely used spare parts based on the BP neural network. Applied Mechanics and Materials, 519520:1513–1519, 2014.
[3] Sadeghi J., Sadeghi S., Niaki S.T.A. (2014). A hybrid vendor managed inventory and redundancy allocation optimization problem in supply chain management: An NSGAII with tuned parameters. Computers & Operations Research, 41:53–64, 2014.
[4] Qiu Z., Ni Z. Probabilistic interval approach for determining the demand of aviation spares. Acta Aeronautica Et Astronautica Sinica, 30(5):861866, 2004.
[5] Axsäter S. Optimization of orderuptos policies in twoechelon inventory systems with periodic review. Naval Research Logistics, 40(2):245–253, 1993.
[6] Kouki C., Babai M.Z., Jemai Z., Minner S. Solution procedures for lost sales basestock inventory systems with compound poisson demand. International Journal of Production Economics, S0925527318300513, 2018.
[7] AlRifai M.H., Rossetti M.D. An efficient heuristic optimization algorithm for a twoechelon (R, Q) inventory system. International Journal of Production Economics, 109(12):195–213, 2007.
[8] Hnaien F., Delorme X., Dolgui A. Multiobjective optimization for inventory control in twolevel assembly systems under uncertainty of lead times. Computers & Operations Research, 37(11):1835–1843, 2010.
[9] Gicquel C., Cheng J. A joint chanceconstrained programming approach for the singleitem capacitated lotsizing problem with stochastic demand. Annals of Operations Research, 264(12):123155, 2018.
[10] Talaei M., Farhang M.B., Pishvaee M.S., et.al. A robust fuzzy optimization model for carbonefficient closedloop supply chain network design problem: a numerical illustration in electronics industry. Journal of Cleaner Production, 113:662–673, 2016.
[11] Shang C., You F. Distributionally robust optimization for planning and scheduling under uncertainty. Computers & Chemical Engineering, 110(FEB.2):5368, 2018.
[12] Shi Y., Boudouh T., Grunder O. A fuzzy chanceconstraint programming model for a home health care routing problem with fuzzy demand. International Conference on Operations Research & Enterprise Systems, 2017.
[13] Rubinstein R.Y, Kroese D.P. Simulation and the Monte Carlo method by Reuven Y. Rubinstein. Journal of the American Statistical Association, 78(382):511512 1983.
Published on 12/01/21
Accepted on 12/11/20
Submitted on 13/07/20
Volume 37, Issue 1, 2021
DOI: 10.23967/j.rimni.2020.11.002
Licence: CC BYNCSA license
Are you one of the authors of this document?