A hybrid meta-heuristics approach for supplier selection and order allocation problem for supplying risks of recyclable raw materials

Article history: Received May 3 2020 Received in Revised Format May 21 2020 Accepted November 29 2020 Available online December, 3 2020 In this work, a mixed-integer linear programming model is formulated to allocate the appropriate orders to the right suppliers for recyclable raw materials. We modify the previous model for the supplier selection and order allocation problem for stochastic demand to cope with the supply risks of recyclable raw materials such as insufficient supply quantity, defective rate, and late delivery. The optimal solution of the mathematical model is the benchmark for small-sized problems. Then, a hybrid meta-heuristic of Particles Swarm Optimization and Grey Wolf Optimization (PSO-GWO) is proposed to search for the best solution for large-sized problems. A real-life case study of a steel manufacturer with two factories in Vietnam is presented to validate the proposed approach. Some experiments have been tested to confirm the performance of the hybrid PSO-GWO approach. © 2021 by the authors; licensee Growing Science, Canada


Introduction
Industrial wastes are recently encouraged to be used in production processes in many countries. Typically, Gharfalkar et al. (2015) introduced a so-called waste utilization priority triangle (Fig. 1). The activities of waste utilization start from the replacement of non-materials to reuse recyclable resources and end up with returns. In industry, the most common way to use wastes in production processes is to recyclable raw materials. These recyclable raw materials are usually outsourced from small collecting companies. However, the supply from these collecting companies is naturally uncertain due to many risks such as insufficient quantity, unsatisfactory quality, and late delivery. In the literature, two approaches are focused on dealing with these factors. The first direction tries to increase the efficiency of collection processes. Many solutions are proposed for waste treatment, disposal, and collection processes to minimize costs and environmental impacts. Some examples are the works of Badran and El-Haggar, 2006;Chang et al., 2012;Anghinolfi et al., 2013;Zhang et al., 2014;Xu et al., 2017. In this study, we focus on investigating the second approach that integrates the waste disposal and reuse processes with the production system to reduce supply risks. Among the most critical issues, the supply risks integrated with the supplier selection and order allocation (SSOA) problem have been intensively investigated in recent research works (Kim et al., 2011;Al-E-Hashem et al., 2013;Kumar et al., 2016;Hamdan and Cheaitou, 2017;Moheb-Alizadeh and Handfield, 2017;Cheraghalipour and Farsad, 2018;etc.). Besides dealing with the uncertain of the demand, different aspects of supply risks have been considered but the most popular ones could be the defective rate, available quantity, and late delivery. In an early work of Kim et al. (2011), some risk factors such as shipping errors, missing, and wrong goods are taken into account. The appropriate order quantity for the retailer has been determined subject to the probability of defective goods. The ordering plan is also frequently updated to cope with uncertain demand.  Gharfalkar et al.,2015) For the varying supplier availability, Hamdan and Cheaitou (2017) solved the multi-period green supplier selection and order allocation problem with quantity discounts and availability of suppliers from period to period. They used a fuzzy TOPSIS and AHP based approach to select the suppliers to input into a single product bi-objective integer linear programming model with deterministic demand. However, only the supplier availability instead of supply quantity is addressed in this study.
Regarding delivery risks, Al-E- Hashem et al. (2013) presented a random planning model to determine supplier orders for a green supply chain integrated with production planning with demand uncertainty. In this model, the correlation between waiting time, shipping costs, waiting time, and greenhouse gas emissions is also considered. Among the popular supply risks of quality, delivery, and supply quantity, a lot of works dealt with the similar risks such as late delivery and disruptions (Sawik, 2014;Tsai, 2015;Moghaddam, 2015;Venkatesan and Goh, 2016;Vahidi et al, 2018;Cheraghalipour and Farsad, 2018). Sawik (2014) presented a stochastic mixed-integer programming model for supplier selection, order allocation, and order scheduling subjects to disruption risks. Both risk-neutral and risk-averse solutions are proposed to determine the expected production schedule. Later, Tsai (2015) and Moghaddam (2015) addressed both rejections/defective parts and late delivery. While Tsai (2015) developed only a mixed-integer non-linear programming formulation, Moghaddam (2015) used a hybrid Monte Carlo simulation integrated with goal programming method to solve the multi-objective model for the supplier selection and order allocation problem in closed-loop supply chain systems. In another work, Venkatesan and Goh (2016) formulated a multi-objective MILP model for selecting the supplier and allocating order quantity under disruption risks. They first evaluated and ranked suppliers by the fuzzy AHP-PROMETHEE approach. Then, the multi-objective PSO algorithm searched the compromised solutions to the problem. Recently, Cheraghalipour and Farsad (2018) also studied the sustainable supplier selection and order allocation problem in multi-period, multi-item, and multi-supplier environment considering quantity discounts and disruption risks. They implemented the so-called best worst method to derive the weight of criteria and suppliers. Then, they revised the multi-choice goal programming to solve the bi-objective MILP model for calculating the amount of inventory and shortage, ordered quantity for each selected supplier. The effects of quantity discounts and disruption risks are also analyzed. Besides disruption risks and uncertain demand, Vahidi et al. (2018) also concerned the operational risks for the sustainable supplier selection and order allocation problem. Some strategies are used together to select a suitable supply base such as a hybrid SWOT-QFD framework, contracting with backup suppliers, and considering excess capacity for supply nodes.
From the above review, few research works investigated the most popular supply risks of defective rate, late delivery, and varying supply quantity together with uncertain demand. Only three works of Govindan and Sivakumar (2016), Kumar et al. (2016), and Moheb-Alizadeh and Handfield (2017) so far addressed these aspects at a certain level. Gonvindan and Sivakumar (2016) proposed a two-phase hybrid approach. The first phase presents the rating and selection of potential suppliers by considering economics (cost), operational factors (quality and delivery), and environmental criteria (recycle capability and GHG emission control) by Fuzzy TOPSIS methodology. The second phase presents the order allocation process using multiobjective linear programming to minimize cost, material rejection, late delivery, recycle waste, and CO2 emissions in the production process. Kumar et al. (2016) represented how to optimize orders among various suppliers while taking into consideration all three dimensions of sustainability, economic, social, and environmental. An integrated fuzzy AHP and fuzzy multi-objective linear programming approach are used to allocate orders among suppliers. Various factors such as quality, lead time, cost, energy use, waste minimization, emission, and social contribution are weighted in the AHP framework. Then, these weights are taken into account in multi-objective linear programming to determine the optimal order allocation with WASTE NON-WASTE fuzzy demand. The recent research of Moheb-Alizadeh and Handfield (2017) formulated a multi-objective mixed-integer nonlinear programming (MINLP) model for the SSOA problem with uncertain demand, deterministic defective date, late delivery quantity, and limit supply quantity.
In this work, we modify the model of Moheb-Alizadeh and Handfield (2017) for stochastic demand, defective rate, late delivery quantity, and varying supply quantity for a case of medium-sized enterprises producing steels from recyclable raw materials used in construction applications. The main difficulties in supplying these recyclable raw materials are the instability of supply capacity, quality, and delivery time. These supplying limitations and unstable demand push up the purchasing cost of the enterprise. Therefore, it is required to design an effective sourcing strategy to reduce the fluctuations of supplying capacity, compensate for the loss of quality and late delivery. Also, an appropriate order volume not only mitigates the supplying risks but also reduces efficiently total operating costs in terms of the purchasing price, shipping costs, and inventory level. To overcome these challenging tasks, a mixed-integer linear programming model for supplier selection and order allocation recyclable raw materials is figured out in the next section. The optimal solution obtained from this model for smallsized problems is used as a benchmark for our meta-heuristic solution developed in Section 3. In Section 4, we investigate a steel manufacturing company in Vietnam to show our efficient results. Some analyses are also conducted in this section to confirm the advantages of our proposed approach. Finally, conclusions and some recommendations for future research are presented in Section 5.

Supplier Selction and Order Allocation Model for Recycle Raw Materials
In this section, we modify the similar model of Moheb-Alizadeh and Handfield (2017) for our supplier selection and order allocation problem for recyclable raw materials. Moheb-Alizadeh and Handfield (2017) introduced a multi-objective mixedinteger non-linear programming model for supplier selection and order allocation problem with stochastic demand. They consider multiple objectives of cost, environmental, and social impact. Here, we modify the model of Moheb-Alizadeh and Handfield (2017) with only cost objective and updating some constraints related to supply risks for recyclable raw materials as follows.

Indices:
: the index of suppliers; = 1, … , : unit purchasing price of raw material j at price level s provided by supplier i in period t.
: fixed shipment cost of raw material j shipped from supplier i to factory k in period t.
: variable shipment cost of unit raw material j shipped from the supplier i to factory k in period t.
ℎ : holding cost of each unit of raw material j at factory k in period t.
: shortage cost of each unit of raw material j at factory k in period t.
: rejection percentage of raw material j provided by supplier i in period t.
: rejecting rate limit of raw material j at factory k in period t.
: on-time delivery percentage of raw material j provided by supplier i in period t : minimum acceptable percentage for on-time delivery of raw material j to factory k in period t.
: required demand of raw material j at factory k in period t.
: maximum quantity of raw material j that supplier i can provide in period t.
: required space for holding a unit of raw material j : maximum available storage space for raw material j at factory k in period t.
: maximum amount of shortage for raw material j at factory k in period t.
: maximum perished amount of material j at factory k in period t : lower limit of purchase volume for raw material i in supplier j at price level s in period t 180 : upper limit of purchase volume for raw material i in supplier j at price level s in period t ̅ : average capacity of transportation vehicles Decision variables: : order quantity of raw material j at price level s provided by supplier i for factory k in period t : inventory level of material j at the end of period t at factory k : shortage level of material j at the end of period t at factory k =1 if raw material j is stored by factory k in period t; 0, otherwise. = 1 if raw material j is lacked by factory k in period t; 0, otherwise. =1 if supplier i provides raw material j at price level s to factory k in period t; 0, otherwise.
The Model: In this model, we consider K factories to buy J (recyclable) raw materials from I suppliers. From now on, we use raw material(s) in short for recyclable raw material(s) with the same meaning. Each raw material will be bought at a certain price level s from its associated supplier i. We consider T periods of the planning horizon. We also assume that the demand, maximum supply quantity of each supplier randomly follows the normal distribution. Purchasing cost includes material cost and shipping cost from the supplier to the factory. The objective function (1) minimizes the total cost including purchasing cost, fixed and variable transportation costs, holding cost, and shortage cost. Constraint (2) enforces that the total rejected amount of each raw material at each factory should not more than the allowable rejection quantity. Constraint (3) ensures that the total on-time delivery quantity of raw materials must meet the required demand at each factory for each period. Constraint (4) is supplying capacity limits. Constraint (5) is an inventory balance equation at each factory. Constraint (6) determines the shortage amount of each raw material at each factory if any. In the beginning, we assume that = = 0. Constraint (7) guarantees the available storage space at each factory for each raw material in each planning period. Constraints (8) and (9) require that the shortage or storage amount of each raw material must not exceed the allowable shortage or a storage quantity, respectively, in each period at each factory. Constraint (10) is a mutually exclusive condition for storage and a shortage quantity of each raw material at each factory in each planning period. Constraints (11) and (12) determine the allowable range of order quantity at each price level l for each raw material. In practice, the lower bound quantity corresponds to the Minimum Order Quantity (MOQ). The upper bound indicates the budget for buying each raw material at each factory in each period. Constraint (13) is the binary and non-negative conditions of variables. The model of Moheb-Alizadeh and Handfield (2017) handled only stochastic demand by using the following chance constraints and Lemma 1: where is the confidence level to the stochastic constraint v, x is an n-dimensional decision vector, w is a stochastic vector Lemma 1. (referred by Moheb-Alizadeh and Handfield, 2017) assuming ( , ) = ℎ( ) − , the deterministic equivalent of ( , In the model (1) -(13), we consider not only random demand but also some risk factors such as instability of supplying capacity, defective quantity, and late delivery for recyclable raw materials. All constraints included in these stochastic values in our case will be expressed in the form of chance constraints. To convert these constraints into the equivalent deterministic ones, let Because , , , and are stochastics parameters, and can be assumed to follow normal distribution with: From the results of chance constraints, we have: From the above results, the following equivalent deterministic MILP model is derived for our supplier selection and order allocation problem subjects to risks of recyclable raw materials: Minimize (1) Subject to: (19) -(23), (7) - (13). This deterministic mixed-integer linear programming model can be solved optimally for small-sized problems by standard software such as LINGO or CPLEX. When the problem sizes become larger, meta-heuristic techniques such as Genetic Algorithm (GA), Simulated Annealing (SA), Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO), Grey Wolf Optimization (GWO) are popular approaches to search for the best solution in a reasonable time.

Solution Approach
As we have seen in the previous section, meta-heuristic approaches have been applied successfully to search near-optimal solutions in a reasonable time for discrete optimization problems. However, a drawback is a difficulty in the selection of the appropriate methodology and its associated parameters for each application. This can be overcome by large-scale numerical experiments. Therefore, we will investigate in this study the popular GWO and PSO techniques that recently show efficient in solving similar applications like our supplier and order allocation problem.

Grey Wolf Optimization (GWO)
Grey Wolf Optimization (GWO) was first introduced by Mirjalili et al. (2014). It imitates the behavior of the wolves in catching their prey. In the GWO, the wolves in the pack are dominated by three key wolves: alpha (α), beta (β), and delta (δ).
The alpha is the fittest solution that takes the role of a leader to make decisions about hunting. The second and third best solutions are beta and delta which play the role of the advisor to support the alpha in making decisions and the role of subordinate to submit to alphas and betas, respectively. The normal wolves of the pack are omegas (ω) which are the rest of the candidate solutions. In the GWO algorithm, the hunting (optimization) process is guided by α, β, δ. The ω wolves follow these three wolves. The hunting process of the wolves has been mathematically modeled the searching, encircling, hunting, and attacking prey (see Fig. 2). The search process starts with creating a random population of grey wolves (candidate solutions). In each iteration, the grey wolves encircle prey by estimating the distance from each candidate solution (wolf) to the prey.

Fig. 2. The Typical GWO Algorithm
where indicates the current iteration X ⃗ is the position vector of the prey X ⃗ is the position vector of a grey wolf.
A ⃗ and C ⃗ are coefficient vectors, in which: Vector ⃗ are linearly decreased from 2 to 0, Vectors ⃗ and ⃗ are random in [0, 1].
The hunt is usually guided by the alpha. The first three best solutions (α, β, δ) obtained so far will estimate the position of the prey and oblige the other wolves to update their positions randomly around the prey (Figure 3).
The parameter a is decreased from 2 to 0 to emphasize exploration (searching for prey) and exploitation (attacking prey), respectively. Candidate solutions tend to diverge from the prey when ⃗ > 1 and converge towards the prey when ⃗ < 1.
In addition, the ⃗ vector provides randomly weights for prey that could help GWO to be more random behavior, favoring exploration, and local optima avoidance.

Particle Swarm Optimization (PSO)
Particle Swarm Optimization (PSO) was proposed by Kennedy and Eberhart (1995). It models the food searching process of the flocks of birds. In the PSO, each particle representing a solution has a position and a velocity at each iteration. The main purpose of the PSO algorithm is to generate a number of particles moving toward the optimal value at a randomized velocity. The PSO keeps track at each iteration the best value of a particular particle (pbest) and the best solution in the neighborhood search (gbest) to guide the search process toward the global optimal solution. The typical PSO pseudo-code is illustrated as follows:

Procedure PSO Begin
Step 1. Initialize parameters of the PSO problem

Do
Step 2. Evaluate fitness value (objective function value) Step 3. Update personal best (pbest) and global best (gbest) values Step 4. Update the inertia weight, cognitive coefficient, social coefficient, the velocity, and the position of each particle

While (termination condition is not satisfied) End
In the PSO, the velocity and position are normally updated by the following formulae: where t ik v : Velocity of particle (i) at time (t) on dimension (k) t ik x : Position of particle (i) at time (t) on dimension (k) ik t pbest x : The particle (i)'s individual best solution as of time (t) on dimension (k) k t gbest x : The Swarm's best solution as of time (t) on dimension (k) 1 c : Cognitive coefficient, it acts as the particle's memory, causing it to return to its individual best region of the search space.
This coefficient limits the size of the step the particle takes toward its personal best.
2 c : Social coefficient, it causes the particle to move to the best regions the swarm has found so far. The coefficient limits the size of the step the particle takes toward the global best.By empirical studies, we often set = = 1 or = = 2.
Here, the position represents the solution. The velocity modifies solutions that help to search for better solutions in the next iteration. The velocity is originally influenced by two coefficients c1 and c2 that determine whether a particle moves toward the best so far position or the global best position of the whole swarm, respectively. However, the values of velocity may always increase in this updating scheme. Shi and Eberhart (1998) introduced an additional parameter of inertia weight w t to control the value of the velocity at the next move (t+1) from the previous iteration t. The inertia weight w t keeps the particle moving in the same direction it was originally heading. Lower values of w t speed up the convergence to the global best, while the higher ones encourage exploring the search space. The complete updating scheme for the velocity will be: wmin and wmax are minimum and maximum inertia weight values, respectively.

Solution Representation
In this section, we present the solution structure of the problem. First of all, we encode the order quantity of raw material j at price level s provided by supplier i for factory k in period t, , as a multi-dimensional non-negative vector ⃗ . In detail, consider a solution of the problem with I = 3 suppliers; J = 2 types of raw material; K = 3 factories to buy recycle raw materials at S = 2 price levels in T = 2 periods. The complete solution vector X will be: Fig. 4

. Solution Representation
In this representation, Xtijsk represents for the order quantity at time t by factory k from supplier i for raw material j at price level s. For example, X23213 indicates the order quantity at period 2 from supplier 3 for raw material 2 at price level 1 for factory 3. This vector X is treated as a particle in our proposed PSO.

Fitness Function
The fitness value is evaluated from the total cost. From (1), we have the total cost including purchasing cost, fixed shipping cost, variable shipping cost, inventory cost, and shortage cost. The steps to calculate the fitness function are as follows: Step 1: Convert solution ⃗ from vector form to multidimensional variable form Step 2: Apply equation (22)

and (23) to calculate inventory and shortage
Step 3: Checking the feasibility of solution ⃗ by capacity and limitation constraints Step 4: If solution ⃗ is not feasible. Update solution ⃗ and Go back to step 2. Otherwise, continue to step 5 Step 5: Use equation (1) to calculate the total cost.

The Proposed Hybrid Meta-heuristics
As we have seen above, the GWO has been controlled by two vectors ⃗ and ⃗ for diverging and converging towards the prey, respectively. In the PSO, the direction of particle flies is influenced by three coefficients of inertia weight w t , acceleration weights c1 and c2. The purpose of setting these parameters is to balance between exploration and exploitation in the searching process. However, these factors are biasedly input by predetermined values. This is not good to help the wolves or particles to adjust to proper searching direction during the iterative process. The wolves or particles have a high tendency to converge fast and premature. This leads to an increase in the chance of falling into local traps.
Besides, it can also recognize that the GWO is strong in converging towards the prey by adjusting the only vector ⃗ < 1. The PSO is easily diversified by changing the acceleration coefficients c1 and c2 in each iteration instead of vectors ⃗ or ⃗ . Moreover, vector ⃗ < 1 serves for exploitation purposes only. Therefore, it is better to combine these two meta-heuristics into one searching procedure to enhance converging speed but still to be able to avoid local traps. In this study, we will update simultaneously the cognitive coefficient , social coefficient , and inertia weight (w t ) of the PSO in each iteration to diversify the searching space. The GWO is just served for converging purposes only if it is prone to stagnation in local situations. The basic steps of the combined meta-heuristics PSO-GWO are outlined in the following pseudo-code.

Procedure PSO-GWO Begin
Initialize the particle population Intitialize parameters while (t< Max number of iterations) for each particle with position xp do a d j u s t xp calculate inventory and shortage quantity check feasibility of xp, inventory, shortage quantity whilexp, inventory, shortage quantity is not feasible calculate fitness value f(xp) iff(xp) is better than pbestp then pbest p ←x p endif iff(pbestp) is better than f(gbest)then gbest←pbestp endif end for update inertia weight w t for each particle with position xp update , calculate velocity of each particle update position of each particle end for if rand (0,1) <prob Initialize the grey wolf population Initialize parameter a, A and C for each search agent (SAg) do adjust (SAg) calculate inventory and shortage quantity check feasibility SAg, inventory, shortage quantity while SAg, inventory, shortage quantity is not feasible Calculate the fitness of each search agent endfor = the best search agent = the second best search agent = the third best search agent while (u< Max number of iterations) for each search agent Update the position of the current search agent endfor update a, A and C do adjust (SAg) calculate inventory and shortage quantity check feasibility SAg, inventory, shortage quantity while SAg, inventory, shortage quantity is not feasible calculate the fitness of each search agent update , , u = u + 1 endwhile xp = mean positions of three search agent , , endif t=t+1 endwhile return gbest End.
In this proposed PSO -GWO algorithm, the quality of solution and computation time are influenced by important parameters such as the acceleration coefficients c1 and c2, the inertia weight w t , the number of particles, the number of iterations of the background PSO, the number of wolves, and the number of iterations of the GWO algorithm. The parameters c1 and c2 represent the learning factor of each particle to diversify the searching space. The inertia weight coefficient w t maintains the particles moving in a certain direction. On one hand, the smaller the inertia weight w t is, the faster the convergence process will be. On the other hand, the greater the value of w t is, the larger the searching space will be. This leads to a larger number of iterations to search for the global optimal solution and longer computation time. Therefore, the GWO algorithm combined with the PSO could help to enhance the search to produce good results. The GWO coefficients are often set at low levels to accelerate the converging process.
Here, the PSO's weights , , and w t are adjusted automatically based on the current status of solutions at each iteration (the global best and local best values): where: ( ) is the fitness of particle at iteration The detailed experiments in the following section will select specifically these parameters. With all settings for the proposed PSO-GWA algorithm, the particles of the PSO will be able to explore the search space to investigate potential good solutions. Then, the GWO could help to converge quickly to the global optimal solution.

Experimental Design
In this section, we will investigate the recyclable raw materials supplier selection and order allocation problem for the case of a steel manufacturer in Vietnam. The company has two factories in Dong Nai province and Da Nang city, Vietnam. The recycle steel raw materials are ordered from collecting units in neighboring provinces. Historical data of demand collected from two factories could approximately follow the normal distribution (Table 1). Due to uncertain market, the suppliers can only commit a certain level of defective rate and late delivery, and quantity of recyclable steel raw materials. Because data has been collected for a long time, other parameters also vary in the form of uniform distribution in our testings as the following table. Confident level for all test cases is 95%. For small-sized problems, we are using the optimal solution of the MIP model as the benchmark to compare with the proposed hybrid PSO-GWO algorithm. The parameters of PSO and GWO are initially selected as follows: wmin = 1.2, wmax = 1.6,c1 = 1.2,c2 = 0.5, number of particles = 20, number of main iterations = 50, number of wolves in the pack = 20, number of GWO's iterations = 10. The sensitivity analysis could be conducted to identify the best combination of these parameters. With the current selection, the experimental tests are conducted on the Xeon E7 CPU and 32GB RAM computer. The results of 10 small-sized test cases are given in the following Table 3. The computation time of the proposed algorithm is about 2000 seconds.  For large-sized problems, we conduct 24 test runs to compare our hybrid PSO-GWO with the stand-alone PSO and GWO for different scenarios by changing the number of suppliers, numbers of raw materials, planning horizon, etc. The PSO, GWO, and PSO-GWO are using the same input parameters. The results are compared in terms of computation time and objective value (Table 4).

Table 4
Experimental Results for Large Sized Problems The obtained results are quite promising. Our PSO-GWO solution is very close to the optimal solution at an average level of 4.5% for small-sized problems. For the large-sized problems, the proposed hybrid PSO-GWO gives better solutions than the stand-alone PSO or GWO, respectively. On average, the gap between GWO and PSO-GWO is 10.56%. The gap between PSO and PSO-GWO is 7.87%. The reason that PSO could give a better solution than GWO comes from the ability to escape from local optima easier by the acceleration weights c1 and c2. The hybrid PSO-GWO is better for both stand-alone PSO and GWO because it utilizes the advantage of each one. Next, we will investigate the effect of rejection rate, on-time delivery rate, supply capacity, and demand on the total cost. Fig. 4 shows the effect of rejection rate on total cost at three cases of U~(1,4) for case 1, U~(4,9) for case 2, and U~(9,13) for case 3. For each case, we conduct the experiments with four confidence levels at 80%, 90%, 95%, and 99%.

Fig. 4. Effect of Rejection Percentage on Total Cost
From Fig. 4, it is observed that the higher the confidence level, the more the total cost will be. With rejection rate at low (case 1), medium (case 2), and high (case 3) levels, when the confidence levels increase from 80% to 99%, the total cost will increase by about 4.04%, 9.98%, and 16.23%, respectively. In addition, the higher the refection rate, the total cost also increases significantly. For example, at the same confident level of 99% the total cost at case 3 (high rejection rate) will be higher than at case 2 (medium rejection rate) about 11.61%, and at case 2 about 22.51 %. These facts come from many reasons. First, if we do not receive the order quantity with the required quality as planned, we need to pay more costs for return, compensation, and additional shipment. In addition, we need to keep more stocks to avoid the higher missing quantity due to a higher defective rate. That explained why the higher probability of a higher defect rate will result in a more increasing percentage of the total cost. Similarly, we conduct experiments to check the effect of the on-time delivery rate on the total cost. In detail, we will change the on-time delivery rate at three levels of U~(92,100) as in practice. The confidence levels are also varying from 80%, 90%, 95%, and 99%.

Fig. 5. Effect of On-time Delivery on Total Cost
From the obtained results in Fig. 5, we can see that the more the certainty of on-time delivery with higher confidence levels, the lesser the total cost will be. When the confidence levels increase from 80% to 99%, the total cost decreases by about 2.94%. These results reflect the fact that if the delivery is as planned (high confident level or high probability of late delivery at the planned rate), we do not need to keep more stock to compensate for the missing quantity due to late delivery. Finally, we are also testing the effect of supply capacity and demand uncertainty on the total cost at four confident levels of 80%, 90%, 95%, and 99%, respectively. Fig. 6 and Fig. 7 show that when the certainty levels of supply capacity and demand increase, the total cost will increase by 18.99% and 11.20%, respectively. It means that if the supply capacity is fluctuated

Effect of On-time Delivery on Total Cost
with high probability due to some reasons such as do not have enough recyclable raw materials to be collected, late raw material supply due to long classification processes or long transportation distance, etc. These changes force the planner to update the plans by keeping more safety stocks, switching suppliers, and placing emergency orders in many suppliers. Moreover, to reduce supply capacity risks, the planner may place orders from many suppliers in different locations which are far away from each other. These could result to increase ordering costs and transportation costs. All these facts lead to an increase in the total cost. Similarly, when the demand is more uncertain, the buyer needs to have some hedging actions to mitigate the risks of demand fluctuation such as a higher level of safety stock, early ordering, etc. In addition, the risk of shortage will be higher. These facts result in higher total purchasing cost. Furthermore, the dependent level of total cost on demand uncertainty is lower than supply capacity (11.2% vs 18.99%) because to mitigate the demand uncertainty, we may increase order quantity until the value of Economic Order Quantity (EOQ). This could increase total purchasing cost but not much as the case of supply capacity uncertainty. In general, the lower the confident levels of rejection rate, supply capacity, and demand; the higher the order volume of each supplier will lead to higher purchasing costs for achieving economic quantity, lower transportation costs because the number of suppliers assigned is low. In contrast, the higher probability of on-time delivery, the lower the total cost because we can avoid keeping more stocks for the missing quantity due to late delivery. Through the above analyses, we could conclude that when the confidence level increases, the total cost also increases due to the insufficient quantity of economic quantity and the increase in the number of assigned suppliers makes the total shipping cost increase.

Conclusions
In this work, we have developed a new mixed-integer model for the supplier selection and order allocation for recyclable raw materials under uncertain demand, varying supply quantity, stochastic defective rate, and late delivery. A real case of steel manufacturer in Vietnam is investigated for the validation of the model. We also develop a hybrid PSO-GWO algorithm to solve the SSOA problem subject all these supply risks together. We design numerical experiments to validate our proposed PSO-GWO algorithm. The obtained results show that our solution is very promising. It is very closed to the optimal solution at an average level of 4.5% for small-sized problems and could give better solutions than standard meta-heuristics at about 10%. We also investigate the effect of different stochastic factors on the performance of the order plan. The obtained results show that the lower the confidence levels of rejection rate, supply capacity, and demand; the higher the total purchasing cost could be. In contrast, the higher probability of on-time delivery, the lower the holding cost of safety stocks for the missing quantity due to late delivery will be. These findings are significant guides for the planner to deal with the risks of recyclable raw materials supply. In the future, we may consider more risks or take more objectives to be included in our model for wider applications.