## 2 Problem Formulation
Therefore, the decision maker aims to maximize \(Z_{1}\) and \(Z_{3}\) and minimize \(Z_{2}\). Equations (7) or (8) ensures that the total transported quantity from \(s\) for \(i\) at \(t\) will be less than, or equal to, on-hand quantity. Equations (9)-(11) are modeled for selecting the k from s to p and not to exceed its capacity where M is a big number. Equation (12) is a stock balance constraint for \(i\) in \(p\) at \(t\). Equation (13) is a capacity constraint for inventory level of \(i\) in \(p\). Equations (14)-(18) are capacity constraints for regular production, overtime production and product inventory levels in p, respectively. Equation (19) is an inventory balance constraint for \(j\) in \(p\). Equations (20)-(22) are constructed to select the \(k\) from \(p\) to the \(w\) and not to exceed its capacity. Equation (23) is an inventory balance constraint in \(w\). Equation (24) is an inventory capacity constraint for a \(w\). Equations (25)-(27) are constructed to select \(k\) from \(w\) to \(r\) and not to exceed its capacity. Equation (28) is a balance constraint for inventory and backorderlevel in \(r\). Equation (29) is an inventory capacity constraint for \(r\). Equations (30) or (31) ensures to meet customer demand at a given probability level. Equation (32) is the definition of the decisions variables. In order to solve PDP, the MOMILP is converted into an equivalent single-objective 0-1 mixed integer programming model by using [PERSON]'s [52] fuzzy programming method. According to the fuzzy programming method, first, the positive ideal solutions (PIS) and negative ideal solutions (NIS) of the objective functions can be specified as Equation (33) [53]:  
\[\begin{array}{l}Z_{1}^{PIS}=min(C^{m}-C^{p})x,\qquad Z_{1}^{NIS}=\max(C^{m}-C^  
For \(Z_{1}^{\text{\emph{NIS}}}\): Max \(Z_{1}\)  
\(s.t.\) Equations (8), (10-(13), (15), (17)-(19), (21)-(24), (26)-(29), (31) and (32)  
For \(Z_{2}^{\text{\emph{PIS}}}\): Max \(Z_{2}\)  
\(s.t.\) Equations (8), (10-(13), (15), (17)-(19), (21)-(24), (26)-(29), (31) and (32)  
For \(Z_{2}^{\text{\emph{NIS}}}\): Min \(Z_{2}\)  
\(s.t.\) Equations (7), (9), (11)-(14), (16), (18)-(20), (22)-(25), (27)-(30) and (32)  
For \(Z_{3}^{\text{\emph{PIS}}}\): Max \(Z_{3}\)  
\(s.t.\) Equations (8), (10)-(13), (15), (17)-(19), (21)-(24), (26)-(29), (31) and (32)  
For \(Z_{3}^{\text{\emph{NIS}}}\): Min \(Z_{3}\)  
\(s.t.\) Equations (7), (9), (11)-(14), (16), (18)-(20), (22)-(25), (27)-(30) and (32)  
_Step 9_: Solve \(Z_{4}\)  
For \(Z_{4}\):  
\(\max\lambda\)  
\(s.t.\)\(f_{i}(z_{i})\geq\lambda,\quad i=1,2,3\) and Equations (8), (10)-(13), (15), (17)-(19), (21)-(24), (26)-(29), (31)  
and (32)  
_Step 10_: If the DM is not satisfied with the initial solution, return to Step 4 and update the \(\alpha\) and \(\beta\) values and repeat the remaining steps until a satisfactory solution is found. The proposed modeling and solution approaches can be solved global optimally by using GAMS's Cplex solver. However, if the numbers of echelons (\(s\), \(p\), \(w\), \(r\), \(k\) and \(t\)) increase, the number of binary variables in the model will increase exponentially. Therefore, optimization packages cannot solve the PDP global optimum, which is a reason to develop a meta-heuristic algorithm.