key: cord-0763923-tbx99q08 authors: Nikzamir, Mohammad; Baradaran, Vahid title: A healthcare logistic network considering stochastic emission of contamination: Bi-objective model and solution algorithm date: 2020-08-25 journal: Transp Res E Logist Transp Rev DOI: 10.1016/j.tre.2020.102060 sha: ef0abccb225bdad3236abc06e3e15983c1a42424 doc_id: 763923 cord_uid: tbx99q08 This paper presents a novel healthcare waste location-routing problem by concentrating on a new perspective in healthcare logistics networks. In this problem, there are healthcare, treatment, and disposal centers. Locations of healthcare centers are known, however, it is required to select appropriate locations for treatment, recycling, and disposal centers. Healthcare wastes are divided into infectious and non-infectious wastes. Since a great volume of healthcare wastes are infectious, the emission of contamination can have obnoxious impacts on the environment. The proposed problem considers a stochastic essence for the emission of contamination which depends on the transferring times. In this respect, transferring times between healthcare and treatment centers have been considered as normal random variables. As transferring time increases, it is more likely for the contamination to spread. Having visited a treatment center, infectious wastes are sterilized and they will no longer be harmful to the environment. This research develops a bi-objective mixed-integer mathematical formulation to tackle this problem. The objectives of this model are minimization of total costs and emission of contamination, simultaneously. Complexity of the proposed problem led the researchers to another contribution. This study also develops a Multi-Objective Water-Flow like Algorithm (MOWFA), which is a meta-heuristic, to solve the problem. This algorithm uses a procedure based on the Analytical Hierarchy Process (AHP) to rank the non-dominated solutions in the archive set. By means of a developed mating operator, the MOWFA utilizes the best ranked solutions of the archive in order to obtain high quality offspring. Two neighborhood operators have been designed for the MOWFA as the local search methods. Extensive computational experiments have been conducted to evaluate the effectiveness of the MOWFA on several test problems compared with other meta-heuristics, namely the Multi-Objective Imperialist Competitive Algorithm (MOICA) and Multi-Objective Simulated Annealing (MOSA). These experiments also include a real healthcare waste logistic network in Iran. The computational experiments demonstrate that our proposed algorithm prevails these algorithms in terms of some well-known performance evaluation measures. healthcare waste management problem. Samanlioglu (2013) presented a multi-objective model for the location-routing problem which has been used for managing industrial hazardous wastes in the Marmara region of Turkey. The proposed model aimed to minimize total costs, transportation risk, and total risk for the population existing around the disposal and treatment sites. Hachicha et al. (2014) modeled the off-site transportation of infectious healthcare wastes as a Capacitated Vehicle Routing Problem (CVRP). Nolz et al. (2014) focused on an inventory routing problem and they applied it for the collection of infectious medical materials. Zhao et al. (2016) developed improved approaches for the network design problem that has been applied in a hazardous waste management system. They have used three multi-objective optimization methods including: (1) Weighted-Sum Method (WSM), (2) Augmented Weighted Tchebycheff Method (AWTM), and (3) Augmented ε-constraint Method (AECM) to solve the proposed problem. Naji-Azimi et al. (2016) focused on the Multi-Trip Vehicle Routing Problem (MT-VRP) for transportation of biomedical samples. They considered a time window for each request and tried to desynchronize the arrivals of biomedical samples through managing: (1) routes ordering, and (2) departure times of vehicles. The researchers developed a multi-start heuristic to solve the model. A grey-AHP approach has been proposed by Thakur and Ramesh (2017) for selecting a healthcare waste disposal strategy in Uttarakhand, Northern State of India. An optimization model has been proposed by Mantzaras and Voudrias (2017) for collection, transportation, treatment, and disposal of healthcare wastes. Hajar et al. (2018) studied the on-site medical waste collection problem and modeled it as a vehicle routing problem with time windows. Rabbani et al. (2018a) proposed a multi-objective model for collection of hazardous wastes. Their proposed model targeted to minimize: (1) transportation and construction costs, (2) population and environmental risks. They developed an augmented -constraint method to solve the problem. In another study, Rabbani et al. (2018b) developed a new model for the industrial hazardous waste locating-routing problem, where there are compatibility restrictions between some types of wastes. To find solutions to their model, they used two meta-heuristics, namely the Non-dominated Sorting Genetic Algorithm II (NSGA-II) and Multi-Objective Particle Swarm Optimization (MOPSO) algorithm. Wichapa and Khokhajaikiat (2018) focused Table 1 Summary of the most related studies. on the location routing problem for the disposal of infectious wastes. They developed a hybrid approach consisting of the AHP, Genetic Algorithm (GA), and Goal Programming (GP) method. Baradaran et al. (2019) studied the stochastic vehicle routing problem, where the fleet of collection vehicles is heterogeneous and there are several prioritized time windows for visiting demand nodes. They developed an Artificial Bee Colony (ABC) method to solve the problem. Babaee Tirkolaee et al. (2019) proposed a mixed-integer linear programming (MILP) model for the multi-trip VRP with time windows in an urban waste collection network. They developed a Simulated Annealing (SA) method to solve the proposed model. Khoukhi et al. (2019) developed a genetic algorithm for the multi-trip VRP with time windows and concurrent pickup and delivery, where a fleet of homogeneous vehicles is used to service a set of hospitals. Ghaderi and Burdett (2019) proposed a two-stage stochastic programming model as an integrated location and routing approach for hazardous materials transportation in multi-modal transportation systems. In the first stage, strategic decisions are made concerning the location of transfer yards and in the second stage, hazmat routing is performed. By reviewing the literature, it has become clear for the authors that the emission of contamination along with other real-life constraints in healthcare logistics transportation can be an interesting topic to focus on. Having reviewed the previous studies, this paper offers the following contributions to the literature: 1. According to the summary of previous studies shown in Table 1 , it is obvious that there is still a need for considering emission of contamination in a healthcare logistic network. Therefore, a mixed-integer mathematical model has been proposed for the healthcare waste location-routing problem considering stochastic transferring times and emission of contamination. 2. To the best of the authors' knowledge, the performance of a water-flow like algorithm has not been investigated on a healthcare waste location-routing problem. Hence, a modified version of this method is developed for the proposed problem. Evaluating the performance of a modified water-flow like algorithm in solving a healthcare transportation problem with the mentioned features was a challenge for the authors that has been tackled in this research. The effectiveness of the proposed algorithm has been demonstrated by means of computational experiments. 3. The proposed algorithm utilizes an MADM-based procedure for generating proper offspring. 4. Two neighborhood operators have been devised for the MOWFA for finding better solutions. This study tries to provide a new mathematical formulation for the healthcare waste location-routing problem. The concern regarding the emission of hazardous contamination during transportation of infectious wastes persuaded us to incorporate this concept into the proposed formulation. The developed model has been designed for two types of wastes which are known as infectious and non-infectious wastes. Both infectious and non-infectious wastes are generated at healthcare centers. Infectious and noninfectious wastes are separated from each other in healthcare centers by trained employees. A heterogeneous fleet of vehicles has been considered to collect infectious and non-infectious wastes, separately. Hence, any interaction between these two types of wastes is avoided. Collection vehicles have limited capacity. The number of available vehicles is predefined and fixed. There are some standards for containers of the vehicles used for collection of infectious wastes. The vehicles for collecting infectious wastes are usually equipped with rigid, tamper-proof, and puncture-resistant containers. These containers are impermeable and made of highdensity plastics or metals that maintain the wastes safely and prevent any liquid leakage (WHO, 2014) . However, these standards are not met in all countries due to high costs and unavailability of the containers with the aforementioned specifications. In the absence of these standards, there is a possibility of emission of contamination to the surroundings. Hence, the proposed model of this paper has incorporated this possibility to minimize the risk of exposure to contamination. In the proposed model, it is not allowed to collect the wastes of a healthcare center partially. The whole amount of each waste type produced at a healthcare center cannot be more than the capacity of the vehicle that is assigned for collection of that type of waste. To collect each type of waste, each healthcare center is visited precisely once by merely one vehicle. The proposed network consists of four main elements: a set of healthcare centers, and three sets of facilities including treatment, recycling, and disposal centers (Rabbani et al., 2018b) . Treatment, recycling, and disposal centers are capacitated. The locations of healthcare centers are known and predefined. However, the locations of treatment, recycling, and disposal centers are determined by the proposed model. A limited number of treatment, recycling, and disposal centers can be established. After visiting a set of healthcare centers, a collection vehicle has two options to drop off its load: (1) if its load is infectious, the vehicle moves to a treatment center, and (2) if its load is non-infectious, the vehicle will visit a recycling center. After treating the infectious wastes, the residues can be either recyclable or non-recyclable. Recyclable wastes are moved to recycling centers, while the others are transferred to disposal centers. Waste residues that are generated by the recycling process will also be transported to disposal centers. The technologies used in treatment centers to treat infectious wastes are not the same. Each treatment center can be established with at most one treatment technology. Time windows are considered for collecting infectious and non-infectious wastes from healthcare centers. Violating a time window will cost a penalty on the network. The time-lapse between generation and treatment of wastes must not exceed 72 h in winter and 48 h in summer. However, due to limitation of space in healthcare centers, both types of wastes must be collected on daily basis. All amounts of infectious and non-infectious wastes must be collected from healthcare centers. All parameters are deterministic except for the traveling times between the healthcare and treatment centers. Graphical display of this network has been depicted in Fig. 1 . The notations used to formulate the proposed model are introduced below: IW : Set of infectious wastes ( = iw I 1, 2, , ). Set of non-infectious wastes ( Recycling cost of the non-infectious waste nw in the recycling centerrc. Based on the assumptions mentioned in Section 3.1, a bi-objective mixed-integer programming model is proposed for the healthcare waste location-routing problem that aims to: (1) determine the routes of vehicles for collecting infectious and noninfectious wastes, (2) Specify the technology that is used by each treatment center, (3) determine the locations of treatment, recycling, and disposal centers, (4) estimate the costs of this network, (5) estimate the amount of emission that will be released to the environment. The formulation of the proposed model is as follows: , , , , , , The second objective function is formulated in Eq. (2) that minimizes the emission of contamination to the environment. Constraint (4) ensures that the amount of wastes (infectious and non-infectious) received by a recycling center cannot be more than its capacity. Constraint (5) guarantees that all amounts of wastes collected by a disposal center cannot exceed its capacity. Constraints (6) and (7) guarantee that all available wastes in this network are collected from healthcare centers. Constraints (8) to (10) indicate that a vehicle must leave the corresponding node of a healthcare center after completion of its service. Constraints (11) to (16) Constraint (17) ensures that the required technology is available in treatment centers for infectious wastes. Equations (18) and (19) guarantee that each healthcare center is visited exactly once. Constraint (20) implies that each treatment center can be established with at most one treatment technology. Constraint (21) to (25) state that the vehicles must be assigned to a facility. Constraints (26) and (27) state that only one vehicle is assigned to the route between healthcare and treatment / recycling centers. Violations of time windows have been obtained in Constraints (28) to (31). Eq. (32) computes the amount of contamination that is released to the environment during transferring the infectious waste iw to treatment centers. Constraint (36) secures that each collection vehicle which has been assigned to carry infectious wastes must drop its load at a treatment center with the compatible treatment technology. Constraints (37) and (38) Equations (39) to (42) Finally, Constraints (43) and (44) (45) and (46), respectively. . The generalized perspective for all infectious wastes proves this lemma. The location-routing problem is known to have an NP-hard nature, which implies that it is complicated to find optimal solutions for large-sized problems (Yildiz et al., 2013) . Besides, the proposed model in this research has a non-linear essence, and therefore it is justifiable to use multi-objective meta-heuristics for solving this problem. Hence, this paper employs two multi-objective metaheuristic algorithms to tackle the problem. In addition to these methods, we have developed a hybrid version of the WFA that differs from the classical version of this method in several perspectives. All employed algorithms are capable of searching for the solution space of the proposed problem in order to find a set of non-dominated solutions for the decision maker to choose from. One of the remarkable factors that influences the complexity and required computational time of evolutionary methods is the representation scheme of solutions (Chen et al., 2013) . Using a proper solution representation may lead the algorithm to all possible candidate solutions (Rabbani et al., 2018b) . To encode the potential solutions of the problem tackled in this paper, we have used the order-based representation scheme which has been used widely in the literature. In this encoding scheme, each solution is a ( matrix, where WT, NF, and NSN represent the number of waste types, centers, and demand nodes, respectively. The first WT rows of a solution are permutations whose elements are in the range ]. The other rows are dedicated to the order of treatment, recycling, and disposal centers. It is possible that these rows have empty cells (Rabbani et al., 2018b) . Suppose that there is a logistic network with dimension (10/5/4/3), which means that this network consists of ten healthcare centers, five treatment centers, four recycling centers, and three disposal centers, respectively. An illustrated example of a solution for this network has been depicted in Fig. 2 . A two-stage decoding procedure is required to transform this representation scheme into meaningful solutions and obtain the objective function values. The first stage is dedicated to the routing dimension of the problem, while the second stage tackles the transportation aspect of the problem. In the first stage, feasible routes for collection vehicles are determined based on their destination nodes. The destination node of a collection vehicle is considered as a treatment, recycling, or disposal center that its remaining capacity is equal to or greater than the capacity of the corresponding vehicle. Having determined the proper center for unloading a vehicle, its collection route will be determined. Hence, an assignment procedure is needed to determine the routes of vehicles. This procedure starts with assigning the first-ranked unassigned node to the route and it continues until no more generation nodes can be added to the route considering the limitations of vehicles. The second stage of the decoding procedure is dedicated to determining the amount of waste residues that should be transferred between different centers. The second stage starts with assigning the residues of the first-ranked treatment center to the first-ranked recycling and disposal centers. This process continues with the rest of the treatment centers according to the ordered list until the remaining capacity of the recycling and disposal centers cannot embrace residues of other treatment centers. The same procedure is applied to transfer non-recyclable wastes from recycling centers to disposal centers (Rabbani et al., 2018b) . Fig. 3 displays the flowchart of the two-stage decoding procedure. Each vehicle passes through these steps regarding the type of its waste. Yang and Wang (2007) save space, the procedures of the classical WFA have been explained in Appendix A. This section introduces the modifications offered in this research for the water-flow like algorithm. One of these modifications has occurred in the splitting and moving procedures of the WFA. In the proposed algorithm, for each original flow, a real random number (RRN) on the interval [0, 1] is generated. Each original flow is substituted by its sub-flows if: (1) the original flow is dominated by at least one of its sub-flows, or (2) the random number RRN is equal to or more than 0.5. Otherwise, the sub-flows are neglected and the original flow is used once again to produce NOSF number of new sub-flows. This procedure forces the MOWFA to find better solutions and enhances the quality and diversity of non-dominated solutions, simultaneously. The first condition helps to detect more appropriate solutions in terms of quality, while the second condition enables the MOWFA to find more diverse solutions. Similar to many other multi-objective optimization methods, the non-dominated solutions found in each iteration of the MOWFA are stored in an archive set. To find the non-dominated solutions of each iteration, the non-dominated sorting technique proposed by Deb et al. (2002) has been used. In this technique, each individual is compared with other individuals of population to determine the first front of nondominated solutions. A multi-objective optimization problem comprises a set of minimizing or maximizing objectives defined over a solution space. For the proposed model with two minimization objectives, a solution j is dominated by another solution j for both objectives and there will be at least one objective such that is less than the one assessed at solution j. The newly found non-dominated solutions are compared to the existing solutions of the archive. A newly found solution becomes a member of the archive if it is not dominated by the solutions of the archive. The dominated solutions of the archive will be eliminated from this set. As another modification to the classical WFA, we propose a new procedure based on the Analytical Hierarchy Process (AHP) and the solutions stored in the archive set. The AHP is a well-known Multi-Attribute Decision Making (MADM) approach proposed by Saaty (1987) to determine the weights of criteria and rank the alternatives. Since the solutions existing in the archive set are the best found solutions during the optimization process, the proposed approach is implemented on these solutions so that more proper solutions might be found in their neighborhoods. In each iteration, the AHP method is applied to the solutions of the archive set, and therefore the weights of the criteria are obtained and the solutions are ranked. The hierarchy used in the MOWFA has three different levels: (1) goal, (2) criteria, and (3) alternatives. Fig. 4 shows an illustration of this hierarchy. The objective functions are considered as the criteria in the AHP method and the solutions of the archive set are considered as alternatives. The relative weights of criteria are determined by means of pairwise comparisons. To conduct these comparisons, the judgments of the decision maker is adopted using the 9-point scale. Keeping the 9-point scale in mind, in the MOWFA, the first objective function has been referred to as "equally to moderately importance" over the second objective function. This means that the first objective function is twice as important as the second objective function for the decision maker. To calculate the weights of criteria, the arithmetic mean as an approximation method is used. The exact calculations have not been considered in the proposed algorithm since the complexity will be increased significantly without any necessity. In the next step, the non-dominated solutions existing in the archive set are compared with each other regarding each objective. The overall priority of alternatives are calculated through combining the weights of objectives with the priority weights of non-dominated solutions. Having ranked the solutions of the archive set, the first and the second ranked solutions are selected for mating. A mating operator is applied to the selected solutions, and therefore two new offspring are generated. It is noteworthy that in each iteration, the Consistency Index (CI) is computed to investigate whether the judgments are consistent or not. If the consistency index is determined to be less than 0.1, the solutions of the archive set are ranked and the best found solutions are selected for mating. Otherwise, the AHP procedure in the MOWFA is neglected for the current iteration and the mating procedure is not performed for the respective iteration. The mating operator generates two binary strings with the length × + WT NF 1 ( ), where WT and NF denote the number of wastes and facilities. Algorithm 1 describes the proposed procedure and the mating operator. Fig. 5 illustrates an example of producing two offspring from the selected solutions. In the next step, two generated offspring are used to update the archive set. ). The nor th element on the first random binary string. The nor th element on the second random binary string. The nor th row of the first selected solution. The nor th row of the second selected solution. The nor th row of the first offspring. The nor th row of the second offspring. Generate two random binary strings; 2. Update the archive set with the newly generated offspring; Eliminate the dominated solutions from the archive set; Output: An updated archive set. In addition to the aforementioned modifications, the MOWFA uses two novel neighborhood operators (as local search methods in the splitting and moving procedures) to produce new feasible solutions from each water-flow. By using these operators, it is possible to find the solutions in the neighborhood of an existing water-flow. These neighborhood operators are elaborated in the following subsections. The MOWFA stops after a maximum number of iterations (MaxIt). Fig. 6 shows the flowchart of the proposed algorithm. A random integer number (RAIN) is generated on the interval [ NSN 1, ] , where NSN denotes the number of demand nodes. The first RAIN columns of a solution is preserved, and several elements of other columns have to take new positions. In this respect, in each row, several elements existing on the columns + RAIN ( 1) to NSN are chosen randomly to change their positions. The pseudocode of this operator is presented in Algorithm 2. First random position number. Second random position number (RNP RNP ) 1 2 . Number of used (non-empty) cells on the nor th row. The value on the RNP 1 column and nor row. The value on the RNP 2 column and nor row. A set to store random numbers. Generate an integer random number ( RAIN 4) . Then, the 5th to 10th columns of the solution are subject to changes in order to produce a new solution. In each row, two positions have been selected randomly and their values have been swapped. To The value on the positionPOS nor 1 . The value on the positionPOS nor 2 . Number of rows in a solution ( = nor NOR 1, 2, , ). A set to store random numbers for the nor th row. Number of used (non-empty) cells on the nor th row. . The values on these positions are swapped, and therefore a new solution is generated as shown in Fig. 8 . To investigate the success of the MOWFA in solving the proposed model, we have used two other meta-heuristics namely, the M. Nikzamir and V. Baradaran Transportation Research Part E 142 (2020) 102060 Multi-Objective Imperialist Competitive Algorithm (MOICA), and Multi-Objective Simulated Annealing (MOSA). The reason behind choosing these two well-known methods is that the MOICA and MOSA have been successful in solving various multi-objective problems (Hosseini and Al Khaled, 2014; Fathollahi-Fard et al., 2019) The imperialist competitive algorithm (ICA), as an evolutionary algorithm, has been inspired by socio-political evolutions. In the MOICA, each solution is considered as a country. Similar to political backgrounds, countries are either imperialist or colonies. Therefore, for each imperialist, there is a set of colonies. The algorithm searches for the best country during the optimization process. The best country has the best possible combination of socio-political features. These socio-political features include economic policy, culture, language, etc (Zandieh et al., 2017) . To obtain different Pareto fronts, the MOICA hires the non-dominated sorting technique after producing initial population of countries. This algorithm utilizes the crowding distance metric to rank individuals of each front. For the MOICA, there is an archive set that stores non-dominated individuals of the first Pareto front in each iteration. By using the non-dominated sorting technique and crowding distance metric, the best countries (solutions) are chosen from the population and they are regarded as imperialists. The rest of the countries are considered as colonies. The power of each imperialist should be evaluated to determine its share of the colonies. In other words, the more the power of an imperialist, the greater its share of the colonies that can be allocated to it. In each iteration of the MOICA, there is an assimilating movement in which the colonies of an imperialist move toward their respective imperialist. Through this movement, the features of colonies will become similar to the features of their corresponding imperialist. Therefore, the colonies of an imperialist are improved. For the MOICA, there is a revolution operator that enhances the exploration ability of the algorithm. This operator selects a random colony with a probability and revolves its position in the search space based on a random manner. Because of assimilation and revolution procedures, if a colony was found to be better than its imperialist, their positions would be exchanged. The imperialists compete with each other by seizing the weakest colonies of one another. The algorithm stops if a predefined termination criterion is met. In this research, an upper bound has been defined for the number of iterations and the algorithm will stop if it reaches this upper bound (Zandieh et al., 2017; Nabipoor Afruzi et al., 2014) . The roots and basic concepts of this algorithm can be found in Atashpaz-Gargari and Lucas (2007). The simulated annealing (SA) is a renowned method inspired by physical annealing process in the field of metallurgy. This method has been developed based on the statistical mechanics concept. In the annealing process, proper low energy state crystallization is achieved by slow cooling of a metal. On the other hand, fast cooling will result in improper crystallization. Gradual reduction of high temperature to low temperature will lead to a crystalline state. Therefore, each temperature must be maintained for a sufficient amount of time so that the crystal will have enough time to reach the minimum energy state. The SA imitates the crystallization process in order to find optimal or near optimal solution. This research employs the Multi-Objective Simulated Annealing (MOSA) which hires the concepts of dominance and Pareto optimality to converge to the Pareto optimal set. In addition to the convergence behavior, the MOSA utilizes the transition palpability to preserve appropriate spread of the solutions throughout the obtained Pareto front. One of the advantages of the MOSA over evolutionary algorithms is that this method requires no big memory to store the population. Hence, the MOSA can find an approximate Pareto front in a short running time (Zidi et al., 2012) . The MOSA encounters two situations in search of the optimal solution. First, the algorithm finds a better solution than the current one. In this case, the better solution is accepted and the current one is discarded. Secondly, the MOSA finds worse solutions that are accepted according to an acceptance probability (Baños et al., 2013) . The acceptance probability of new solutions when they are worse than the current solution depends on the temperature. The acceptance probability is low if the temperature is low. Conversely, in case of being worse than the current solution, the probability of accepting a new solution is high when the temperature is also high (Yannibelli and Amandi, 2013) . The acceptance probability is determined based on the Boltzmann Probabilistic Distribution (Nino et al., 2012) . This acceptance helps the MOSA to dodge the local optimum. The MOSA stores the non-dominated solutions in an unlimited size archive (Zidi et al., 2012) . The newly found solution can be added to the archive if it is not dominated by any of the existing solutions in the archive set (Baños et al., 2013) . The MOSA continues the optimization process until the temperature cools down below a predefined threshold. This section provides several numerical experiments to evaluate the performances of the algorithms. The location-routing problem (LRP) is extremely complex to solve due to the fact that it belongs to the class of NP-hard problems. This complexity increases when the number of demand nodes or the candidates for location facilities are large. Using meta-heuristics for solving the problem is justifiable since the proposed model is more intricate than the classical LRP. Therefore, meta-heuristics are welcome to quickly find feasible solutions to large size instances of this problem (Prodhon and Prins, 2014). We used the MATLAB software (R2015b) to code the algorithms and we ran these methods on an Inter Corei3 personal computer (3.4 GHz CPU) with 8 GB of memory (RAM). To evaluate the MOWFA and compare its performance with the MOICA and MOSA, 45 test problems with different sizes are produced. The test problems 1 to 15 are small-sized problems, while the test problems 16 to 30 are medium-sized problems. The test problems 31 to 45 are considered as large-sized problems. The main parameters of the proposed problem and their ranges have been shown in Table 2 . The performance of a meta-heuristic is affected by the values of its control parameters (Afshar-Nadjafi et al., 2017) . Hence, in this section, the control parameters of the MOWFA are calibrated via the Taguchi Method (TM) to improve the efficiency of this algorithm. By using the Taguchi method, a large amount of information is obtained by conducting the least number of experiments (Hasçalık and Çaydas, 2008) . In this study, a five-level Taguchi design has been considered for the parameters of the MOWFA. A total objective function (a convex combination of objectives) has been used as a response for the Taguchi Where, OFV j 1 and OFV j 2 denote the values of the first and the second objective functions for the j th solution, respectively. The best acquired values of the first and the second objectives are represented as OFV 1 andOFV 2 , respectively. For the solution j, the total objective function (TFV j ) is computed by Eq. (49). (0 1) represents the relative importance of the objectives. In the next step, the total objective function values are converted into the Relative Percentage Deviation (RPD) which is computed by the Eq. (50) (Hosseinian and Baradaran, 2019a): Where, AL represents the total objective function value obtained by a hired algorithm, and B is the best value of the total objective function. The levels which have been considered for the parameters are shown in Table 3 . The chosen levels of parameters are reported in Table 4 . Besides, Fig. 9 depicts the average S/N ratio plots for the small, medium, and large size test problems. Other algorithms have been tuned by the same procedure. The literature of multi-objective optimization problems provides numerous metrics to evaluate the strengths and weaknesses of multi-objective optimizers. This research hires several comparison metrics to evaluate the performances of the MOWFA, MOICA, and MOSA in obtaining more appropriate solution sets. These metrics are described as follows: We have employed two indexes introduced by Li et al. (2014) to compare the algorithms in terms of computational time. These metrics are denoted as ATO and ATA which represent the average computational time to find one solution, and the average computational time to obtain all solutions, respectively. These criteria are computed using the following formulas: Where, NSOL rn is the number of obtained solutions in the rn th run. RN represents the number of runs and Comt rn denotes the computational time in the rn th run. The error ratio (ER) has been used in many multi-objective optimization studies. The ER is used to measure the non-convergence of an approximation Pareto front found by a multi-objective optimization algorithm towards the Pareto optimal front. This metric is calculated by the following formula (Yazdani et al., 2019) : Where, NODS represents the number of non-dominated solutions on the approximation front obtained by an algorithm. e j is equal to 1 if the j th non-dominated solution found by a given algorithm belongs to the Pareto optimal front. Otherwise, e j is equal to 0. Smaller ER values indicate better performance of a method (Yazdani et al., 2019; Hosseinian, et al., 2019; Veldhuizen, 1999) . The spacing metric (SM), developed by Schott (1995) , is a diversity-based performance measure that provides information on the distribution of solutions existing in an approximation Pareto front (Rabbani et al., 2018b) . In a comparison of two Pareto sets, the approximation front with the lower SM value has a better uniformity of distribution of solutions. The SM is computed by the Eqs. (54) to (56) (Yazdani et al., 2019) : Where, SNDS is the set of non-dominated solutions on the Pareto front. d j denotes the minimum distance between the solution j on the approximation front and other solutions. d is the mean of all d j s. OBJ is the number of objective functions in the optimization problem. OFV obj j represents the value of objective function obj for the solution j. Zitzler (1999) proposed the maximum spread (MS) metric, known as diversification metric, to measure the spread of the approximation Pareto front. The normalized MS metric is obtained by the Eq. (57) (Zitzler, 1999; Helbig and Engelbrecht, 2013; Yazdani et al., Tables 5 to 7 detail the computational results of the MOWFA, MOICA, and MOSA in terms of the ER, SM, and MS metrics. Each algorithm has been run for ten times and the average of the obtained values have been reported in these tables. According to the obtained results, the MOWFA has better convergence than the MOICA and MOSA since the MOWFA acquired fewer ER values. Regarding the uniformity of non-dominated solutions obtained by these three methods, it is inferred that the MOWFA outperformed the other methods with respect to the SM metric. Comparing the MS values indicates that the MOWFA has behaved better in producing more diverse solutions than other methods. The performances of algorithms have been compared in terms of their average computational times in Figs. 10 and 11 . These figures show that the MOSA has the least CPU times in finding one solution (ATO metric) and all solutions (ATA metric). The MOWFA is slower than other methods due to the local search strategies embedded in this algorithm. The performances of the three algorithms have been statistically evaluated as well. The statistical comparisons have been implemented using the 2-sample t-test and the following hypotheses. These comparisons have been made in terms of each performance metric at a 95% confidence level. The null hypothesis (H 0 ) assumes the equality of the means of methods regarding a performance measure which implies that there is no significant difference between them. The alternative hypothesis (H 1 ) presents the opposite assumption. Thus, the null hypothesis (H 0 ) is rejected if the p-value is less than 5% . The means of a performance measure value obtained by the algorithms are not significantly different. The means of a performance measure value obtained by the algorithms are significantly different. Each algorithm has solved the test problems for ten times. The following sub-sections discuss the statistical comparisons between the algorithms in terms of the performance metrics mentioned in Section 5.3. • MOWFA versus MOICA Table 8 compares the performances of the MOWFA and MOICA. Results reported in Table 8 indicate that the p-values are less than 5% for the ER, SM, and MS performance measures. Hence, the null hypotheses (H 0 ) for these metrics are rejected which means that the superiority of the MOWFA is evident over the MOICA. On the other hand, the p-values for the ATO and ATA metrics are not less than 5%. Therefore, it can be concluded that the MOWFA is not significantly slower than the MOICA. The MOWFA and MOSA have been compared statistically, and the results have been summarized in Table 9 . Based on the outputs, the p-values are less than 5% for most of the metrics. Therefore, the null hypotheses for these measures are rejected. This means that the MOWFA significantly overcomes the MOSA regarding the ER, SM, and MS. Moreover, it can be inferred that the MOWFA is significantly slower than the MOSA. In the following, we used a two-phase approach based on the Multi-Attribute Decision Making (MADM) methods to rank the algorithms in terms of their results. The goal of this approach is to find the most efficient algorithm in solving test problems. In this respect, the algorithms and performance evaluation metrics are considered as the alternatives and criteria, respectively. In the first phase of this approach, the weights of the criteria are calculated via the Best-Worst Method (BWM), and the second phase involves prioritizing the algorithms by means of the Simple Additive Weighting (SAW) method. This approach is called the BWM-SAW that has been described in Algorithm 4. Decision matrix (Rows and columns represent the algorithms and criteria, respectively). The element on the na th row and nc th column of the decision matrix. , Normalized element on the na th row and nc th column. Ultimate score of the na th algorithm. Phase 1 (BW-M) 1. Determine the most important (the best) criterion. Determine the least important (the worst) criterion. 3. Use a number on the interval [1, 9] to show the preference of the most important criterion over the other criteria. Hence, the Best-To-Others vector (BTO) is formed as is the preference of the criterion nc over the worst criterion W. 5. Determine the optimal weights of criteria by solving the following model (Rezaei, 2015) : Min CONV (58) Subject to: Output: Ranking of algorithms. Tables 10 and 11 report the results of implementing the BWM-SAW methodology for ranking the algorithms. Based on Table 10 , the ER has been identified as the most important criterion, while the ATO has been discovered as the least important attribute. Table 11 clarifies that the MOWFA prevails other methods in all test problems with different sizes. The MOICA and MOSA are placed in the second, and third positions. The -constraint method (Khodemani-Yazdi et al., 2019) was hired to validate the outputs of the algorithms. This exact algorithm solved the small size problems and its results were compared with the MOWFA, MOICA, and MOSA. Fig. 12 shows that the MOWFA obtains closer results to the -constraint method. Therefore, it is concluded that the MOWFA can obtain near-optimal solutions. Some experiments have been conducted to assess the performance of the MOWFA in comparison with the classical WFA. Since the WFA is a single-objective method, it has been run twice for each of the test problems. In each run, one of the objectives used in the mathematical model has been optimized by the WFA. Fig. 13 shows boxplots comparing the best objective function values obtained by the MOWFA and WFA for small, medium, and large size test problems. According to Fig. 13 , the outcome of the comparisons demonstrates that the MOWFA is more competitive than the classical version of this algorithm. As mentioned in Section 4.3, the MOWFA uses the AHP algorithm in each iteration to find the best solutions of the archive for the mating operator. Since the implementation of the AHP requires pairwise comparisons between non-dominated solutions, it is required to assess the impact of using the AHP on the computational times. Therefore, we have conducted an analysis which evaluates the computation efforts with respect to the number of existing non-dominated solutions in the archive set. Three test instances with different sizes have been selected randomly. From twenty number of runs considered for the MOWFA, we have chosen three runs to investigate the impact of the AHP on the computational time of each iteration. These runs are associated with situations where all judgments are consistent. Figs. 14 to 16 show the computational times of each iteration for the chosen three test problems in three runs. Each figure is associated with one of the problems. As shown in these figures, the required computational times of each iteration increases as the number of non-dominated solutions in the archive grows. It can be concluded that the size of the archive set considerably increases the computational times of the MOWFA. In this section, a real-life logistic network is considered which consists of 15 healthcare centers in Tehran Province, Iran. There are 5, 5, and 4 potential locations for treatment, recycling, and disposal centers, respectively. The locations of healthcare centers are known and predefined. A private transportation company is responsible for collecting wastes of healthcare centers. A team of planners are working in this company that determine the routes of vehicles considering the network's constraints. This company must consider the time windows of these centers since the violation of time windows will result in severe penalties. Through waste collection services of this company, the wastes of all healthcare centers are collected and transported to the required facilities. Nonrecyclable wastes are carried to municipal authorized dump (landfill) sites. There are two types of vehicles for collecting wastes. The company uses rear-loader and manual side-loader vehicles for collecting infectious and non-infectious wastes, respectively. Both types of vehicles need a crew of one or two workers to operate properly. In most cases, one of the workers drives the vehicle, while the other one monitors and fills the garbage bed. The capacity of vehicles for collection of infectious and non-infectious wastes are approximately 2.5-3.5 and 4.0-5.5 tons, respectively. The vehicles initiate the process of visiting healthcare centers at 9:00 pm. As considered in the proposed model, it is possible to establish a treatment center with one of the following sterilization technologies: Baradaran Transportation Research Part E 142 (2020) 102060 (1) Autoclaves: These devices are applied for medical purposes to carry out sterilization through elevated temperature and pressure. (2) Gamma irradiation machines: Sterilization through gamma irradiation is one of the most well-known ways to neutralize contaminations. Gamma irradiation machines wipe out the pathogens via breaking down bacterial deoxyribonucleic acid (DNA) which leads to inhibition of bacterial division. (3) Incineration devices: These machines treat waste materials by combusting organisms to ash. Therefore, all organic elements in biohazard wastes are incinerated and sterilized. In this technique, the volume and weight of wastes are decreased about 90% (Alvim-Ferraz and Afonso, 2003; Manga et al., 2011) . (4) Ethylene oxide equipment: Ethylene oxide is a poisonous gas for treating infectious wastes in low temperatures. Ethylene oxide sterilizers invade the nucleic acids of microorganisms and cellular proteins. These equipment are usually used for heat/moisturesensitive wastes. The costs of applying the aforementioned technologies are different. This network has five types of infectious wastes such as materials used for surgeries and laboratory experiments, clothes and bandages collected from patients, body fluids, etc. Moreover, two non-infectious materials are gathered from healthcare centers such as plasters, blister packs for medicine, expired drugs, etc. These wastes are not hazardous and they do not need any special treatment technology to be disinfected. All types of infectious wastes in this network can be sterilized by the aforementioned treatment technologies. All healthcare centers produce both infectious and non-infectious wastes. Table 12 shows the daily amount of infectious (AIW) and non-infectious (ANW) wastes generated in healthcare centers. To space spaces and avoid repetition of centers' name, an identification (ID) has been assigned to each healthcare center. In this network comprising 15 healthcare centers, total amounts of infectious and non-infectious wastes are approximately 18 and 27 tons, respectively. For the infectious wastes, transferring times between healthcare and treatment centers have been considered as normal random variables (with a mean of 60 min and a standard deviation of 10 min). Through consulting with the experts of the company, the transferring times between healthcare and treatment centers have been extracted and investigated to test whether the acquired data follow the normal distribution or not. To conduct the normality test, the Anderson-Darling statistic has been used at the significance level of 0.05. The null hypothesis of this test assumes that the specified data follow the normal distribution. Therefore, if the P-Value is less than 0.05, it can be inferred that the collected data do not follow the normal distribution. Otherwise, the transferring times have been normally distributed. The normal distribution plot of transferring times between healthcare and treatment centers is shown in Fig. 17 . As shown in Fig. 17 , an approximate straight line has been drawn via the plotted points. Since the = P Value 0.399, the null hypothesis is not rejected, and therefore we can conclude that the transferring times follow the normal distribution. Distances between facilities have been extracted from Google Earth. According to the discussions held by the experts, the emission coefficients of infectious wastes ( ) iw ( = iw 1, , 5) range from 0.30 to 0.35. It is possible to establish three treatment centers, two recycling centers, and two disposal centers. Due to large scope of distances between healthcare centers and other facilities, a part of the real-life logistic network has been shown in Fig. 18 . The proposed model has been applied to this network in order to determine the locations of the treatment, recycling, and disposal centers in addition to determining the order of visiting healthcare centers to collect their wastes. The MOWFA, MOICA, and MOSA Fig. 16 . Impact of the AHP algorithm on computational times of the MOWFA for a large-size problem Table 12 Daily amount of wastes generated in healthcare centers (in tons). have been employed to solve the model considering the conditions of this real-life network. Once again, the performance evaluation metrics have been used to compare these methods. To have a robust comparison, each algorithm was run twenty times and the outputs have been shown in Table 13 . Based on the outputs, it is obvious that the MOWFA has achieved considerably better results regarding the ER, SM, and MS metrics. The MOSA has been the fastest method comparing to other methods. The performances of the employed algorithms have been statistically compared as well. Table 14 shows the statistical analyses for the real-life network. According to the outputs, the null hypothesis is rejected for the ER, SM, and MS metrics which indicate the significant superiority of the MOWFA over the MOICA and MOSA. Fig. 19 illustrates the approximate Pareto fronts obtained by these three methods for one of the runs of problem. Comparing to other methods, the solutions of the MOWFA have better qualities (convergence) which implies that the approximate Pareto front of the proposed algorithm strongly dominates the non-dominated solutions of other algorithms. The solutions of the MOWFA have also covered a wider area of the objective space which shows their higher diversity. According to Fig. 19 , the uniformity of the MOWFA's solutions throughout the Pareto front is far better than the uniformity of the Pareto fronts discovered by the MOICA and MOSA. This was expected due to better values of spacing metric. As mentioned earlier, the problem has been solved for twenty times and the order of visited healthcare centers acquired by the MOWFA for a single run has been reported in Table 15. Table 15 also details the capacity of available vehicles and the weights of wastes they carry to their destination. Based on the information provided in Table 15 , the company requires 8 and 7 vehicles for Where, PRE rn nc represents the percentage relative error between the solutions of the algorithm and experts' estimation in the rn th run for the objective functionnc. ALG nc rn is the value of nc th objective function in the run rn by the optimizer. Exp nc is the estimation value of nc th objective declared by the experts of the company. The values given by the experts of the company have been estimated based on their experiences of similar networks that they have previously encountered. Figs. 20 and 21 depict the gaps for the first and the second objective functions, respectively. According to the PRE values reported in Figs. 20 and 21, it can be inferred that the outputs of the MOWFA is different from the given estimations of experts for less than 10%. This implies the high accuracy of the proposed method and it shows that the MOWFA can offer reliable results for the studied network. Regarding the site selection perspective of the proposed problem, Table 16 summarizes the selected sites for treatment, recycling, and disposal centers. Based on the results in Table 16 , the following outlines have been obtained. By providing these outputs to the experts of the company, they approved the order of visited healthcare centers and objective values. 1. For most of the runs, the combination of treatment site "5″ and technology "1" has been selected as the best choice. The MOWFA has chosen this combination in 70% of runs. The fifth treatment site along with the first technology is the most frequent selected combination by the MOICA as it has been chosen in 45% of runs. For the MOSA, this combination has been selected for 55% of runs. 2. The candidate location "3″ has been chosen as the best site for recycling center selected by all algorithms. This site has been selected as the most appropriate location for more than 60% of runs. This research proposed a bi-objective mixed-integer mathematical model for the healthcare waste location-routing problem. The proposed model concentrates on stochastic emission of contamination during transportation of infectious wastes collected from healthcare centers. This phenomenon that may occur in real-world healthcare logistic networks and other realistic assumptions of the proposed model are the major contributions of this research. The objectives of this model are minimization of total costs and emission of contamination. This study also offers a multi-objective water-flow like algorithm to solve various test instances of this highly intricate problem. These test instances include small, medium, and large size problems. In addition to the proposed algorithm, two other meta-heuristics were hired to solve test problems. The results of all three algorithms were compared based on several performance evaluation metrics. Experimental results demonstrate that the MOWFA is more capable than the MOICA and MOSA in terms of the ER, SM, and MS. The model has been applied to a real healthcare waste logistic network in Iran to show its practicality. All algorithms were used to solve this case study as well. Due to the comparisons made between the algorithms, the MOWFA showed a high level of efficacy in comparison with other state-of-the-art algorithms. The following directions can be considered for future studies: 1. Uncertainty in other parameters of the model can also be considered to present a more realistic problem. 2. The developed local search operators can be embedded in other optimizers and their results should be compared with the outputs of this research. 3. Other performance evaluation metrics such as generalized distance, inverted generalized distance, set coverage, and hyper- volume can be used to assess the efficiency of the proposed method. 4. Comparing to single-solution algorithms, the proposed method may need more computation times. Therefore, developing some procedures that reduce the running time of the MOWFA is an opportunity for further studies. 5. Considering the occurrence of accidents during transferring infectious wastes and the spread of contaminants is an interesting subject for future studies. 6. It is recommended to develop new techniques for expediting the MOWFA in the AHP procedure since the calculation of weights requires more computation efforts as the number of solutions in the archive set increases. Mohammad Nikzamir: Data curation, Writing -original draft, Software, Validation, Visualization, Formal analysis, Resources. Vahid Baradaran: Conceptualization, Methodology, Writing -original draft, Supervision, Visualization, Investigation, Data curation, Project administration, Resources. The classical WFA method uses the following procedures to solve an optimization problem. These procedures continue until a termination condition is reached (Yang and Wang, 2007; Nourmohammadi et al., 2019): • Splitting and moving procedures: As the water-flows pass the terrains, they are divided into multiple sub-flows. Then, they tend to flow through lower terrains. These mechanisms are simulated by the splitting and moving procedures in the WFA. Suppose that there is a flow (F ) that is running with the mass (MS ) F and the velocity (VL ) F . The splitting procedure is executed for the flow F if ( × MS VL ) F F exceeds a specific threshold (TR). The number of flows (NOF) starts growing as the splitting procedure carries on. This growth is constrained by defining an upper limit for the number of sub-flows (NOSF) that originated from each flow. The number of sub-flows derived from the flow F is obtained by Eq. (69) (Yang and Wang, 2007; Nourmohammadi et al., 2019) : Where, NSF F represents the number of sub-flows originated from the F th flow. By using the functionint (), the closest integer number to a value is obtained. Having obtained the number of sub-flows originating from each flow, the positions of the newly generated sub-flows must be determined. This procedure is known as the moving process in the WFA. The moving procedure is implemented by using one or more local search strategies that may lead the algorithm to find more appropriate solutions in the neighborhood. Local search methods are iterated for a certain number of repetitions. It should be noted that when × < MS VL TR ( ) F F , there will be no splitting process and the flow only leads to a feasible neighborhood solution. Having implemented the splitting and moving procedures, the original flow is neglected since it has been substituted by the sub-flows. The mass and the velocity of each sub-flow are calculated by Eqs. (70) and (71), respectively (Yang and Wang, 2007; Nourmohammadi et al., Where, MS F ms , and OFV ms are the mass and fitness value of the flowms. VL F ms , is the velocity of the sub-flow ms streaming from the flow F . denotes the gravitational acceleration. OFV F ms , represents the difference between the fitness values of F and ms. • Merging procedure: There is a possibility that flows revisit each other when they are passing through lower terrains. In this case, these flows will be merged using a merging procedure. Suppose that the flow F encounters the flow ms in the same location which implies that these two flows have the same fitness values. The WFA discards the flow ms and the velocity and mass of this flow are merged with the flow F by using the following equations (Yang and Wang, 2007; Nourmohammadi et al., • Evaporation procedure: The water-flows evaporate gradually as they pass the terrains. This characteristic enables the waterflows to develop their search space in order to avoid the local optima trap. To simulate this procedure, the WFA uses an evaporation operator every ET iterations on all flows. By using Eq. (74), a new mass for the flow F after the evaporation process (MS ) • Precipitation procedure: The WFA applies two sorts of precipitation procedures to simulate the phenomenon of returning the evaporated water to the terrain. These precipitation procedures include regular and irregular precipitations. The WFA implements the regular precipitation along with the evaporation process every ET iterations. In the regular precipitation, the evaporated water will be precipitated and therefore the new flows will have different positions. The mass of the flow ms after the precipitation process (MS ) ms  is obtained as follows (Yang and Wang, 2007; Nourmohammadi et al., • The irregular precipitation occurs when all flows of the same generation have zero velocity. In this case, all flows are evaporated and then the same number of flows will return to the terrains by means of the rainfall. Therefore, the newly generated flows may dodge local traps by deviating from their origins. The masses of these newly generated flows are calculated as follows (Yang and Wang, 2007; Nourmohammadi et Where, NOF denotes the current number of flows. The newly generated flows obtained by the regular and irregular precipitations are added to the population of flows. Hence, a remarkable number of flows may have the same fitness values. In this respect, after implementing each precipitation procedure, the merging process is carried out so as to avoid excessive searches. Project scheduling for minimizing temporary availability cost of rental resources and tardiness penalty of activities Improvement and modification of the routing system for the health-care waste collection and transportation in Istanbul. Waste Manage A new model for the hazardous waste location-routing problem Incineration of different types of medical wastes: emission factors for gaseous emissions Imperialist competitive algorithm: An algorithm for optimization inspired by imperialistic competition Developing an applied algorithm for multi-trip vehicle routing problem with time windows in urban waste collection: A case study A Simulated Annealing-based parallel multi-objective approach to vehicle routing problems with time windows Stochastic vehicle routing problem with heterogeneous vehicles and multiple prioritized time windows: Mathematical modeling and solution approach Effect of solution representations on Tabu search in scheduling applications A Fast and Elitist Multi-objective Genetic Algorithm: NSGA-II Medical waste generation and management in medical clinics in South of Iran A green home health care supply chain: New modified simulated annealing algorithms An integrated location and routing approach for transporting hazardous materials in a bi-modal transportation network Routing system for infectious healthcare-waste transportation in Tunisia: A case study Onsite medical waste multi-objective vehicle routing problem with time windows Optimization of turning parameters for surface roughness and tool life based on the taguchi method Performance measures for dynamic multi-objective optimisation algorithms A survey on the Imperialist Competitive Algorithm metaheuristic: Implementation in engineering domain and directions for future research An Energy-efficient Mathematical Model for the Resource-constrained Project Scheduling Problem: An Evolutionary Algorithm A multi-objective multi-agent optimization algorithm for the community detection problem Modeling of the time-dependent multi-skilled RCPSP considering learning effect: An evolutionary solution approach P-GWO and MOFA: two new algorithms for the MSRCPSP with the deterioration effect and financial constraints (case study of a gas treating company) Solving a new bi-objective hierarchical hub location problem with an M/M/C queuing framework A genetic algorithm for solving a multi-trip vehicle routing problem with time windows and simultaneous pick-up and delivery in a hospital complex Alternatives for treatment and disposal cost reduction of regulated medical wastes A discrete artificial bee colony algorithm for the multi objective flexible job-shop scheduling problem with maintenance activities Healthcare waste management in Cameroon: A case study from the Southwestern Region An optimization model for collection, haul, transfer, treatment and disposal of infectious medical waste: Application to a Greek region. Waste Manage Management of hazardous medical waste in Croatia Sensitivity Analysis of Simple Additive Weighting Method (SAW): The Results of Change in the Weight of One Attribute on the Final Ranking of Alternatives A Multi-Objective Imperialist Competitive Algorithm for solving discrete time, cost and quality trade-off problems with mode-identity and resource-constrained situations Location-routing: issues, models and methods A practical vehicle routing problem with desynchronized arrivals to depot A Novel, Evolutionary, Simulated Annealing inspired Algorithm for the Multi-Objective Optimization of Combinatorial Problems A stochastic inventory routing problem for infectious medical waste collection A water-flow like algorithm for solving U-shaped assembly line balancing problems Improved route planning and scheduling of waste collection and transport A New Multi-objective Mathematical Model for Hazardous Waste Management Considering Social and Environmental Issues Using metaheuristic algorithms to solve a multi-objective industrial hazardous waste location-routing problem considering incompatible waste types Best-worst multi-criteria decision-making method The analytic hierarchy process-what it is and how it is used A multi-objective mathematical model for the industrial hazardous waste location-routing problem Fault tolerant design using single and multi-criteria genetic algorithms optimization Geographic information system-based healthcare waste management planning for treatment site location and optimal transportation routing A routing and scheduling system for infectious waste collection Healthcare waste disposal strategy selection using grey-AHP approach Multiobjective evolutionary algorithms: classifications, analyses, and new innovations. PhD dissertation. Air Force Institute of Technology Solving a multi-objective location routing problem for infectious waste disposal using hybrid goal programming and hybrid genetic algorithm Medical waste management -A review Safe management of wastes from health-care activities Water flow-like algorithm for object grouping problems Hybridizing a multi-objective simulated annealing algorithm with a multi-objective evolutionary algorithm to solve a multi-objective project scheduling problem Evolutionary algorithms for multi-objective dual-resource constrained flexible job-shop scheduling problem Planning for meals-on-wheels: algorithms and application Flexible job shop scheduling under condition-based maintenance: Improved version of imperialist competitive algorithm A hybrid NSGA-II-DEA method for the economic-statistical design of the C-control charts with multiple assignable causes A multi-objective simulated annealing for the multi-criteria dial a ride problem Evolutionary Algorithms for Multi-objective Optimization: Methods and Applications Improved approaches to the network design problem in regional hazardous waste management systems Fig. 19 . Pareto fronts discovered by the algorithms