Adaptive Multi-Objective Algorithm for the Sustainable Electric Vehicle Routing Problem in Medical Waste Management

Medical Waste Vehicle Routing Problem

The scientific study of waste management emerged in the late 1960s ( 18 ). The waste collection vehicle routing problem (WCVRP) was initially defined by Kim et al. ( 19 ), addressing practical aspects like multiple disposal site trips and driver lunch breaks. They proposed an extended version of Solomon’s ( 20 ) insertion algorithm to solve this problem.

The study of MWVRP, a variation of WCVRP, has attracted considerable interest in recent years as a result of COVID-19. Unlike the classic WCVRP, MWVRP requires the evaluation of infection risk. Taslimi et al. ( 21 ) presented a periodic load-dependent VRP for the collection of medical waste, using a decomposition-based algorithm with column generation. Ghannadpour et al. ( 22 ) developed a sustainability-focused health care waste collection model, minimizing transportation costs, fossil fuel consumption, and disease risk, and optimized it with a self-adaptive evolutionary algorithm. Babaee and Aydın ( 23 ) developed a multi-depot MWVRP model considering outsourcing strategy and used meta-goal programming to address multi-objective challenges, including social risk, measured by the population near disposal sites and disease transmission probability. In another study, Tirkolaee et al. ( 24 ) optimized a multi-period MWVRP in uncertain demand situations using fuzzy chance-constrained programming and a weighted goal optimization approach. Govindan et al. ( 25 ) and Zhao et al. ( 26 ) addressed uncertain demand using fuzzy goal programming and augmented ϵ -constraint approaches, respectively. They embedded several practical elements into their models, including facility location planning, vehicle scheduling, and population risk.

In contrast to the aforementioned studies, a recent study by Erdem ( 7 ) incorporated EVs into medical waste collection operations. They focused on reducing health risks and providing cost-effective and zero-emission logistics through a weighting method that combines ALNS and LS as the search engine.

As depicted in Table 1, recent studies on MWVRP have identified various objectives and constraints related to sustainability. The majority of research predominantly focuses on internal combustion engine vehicles (CEV). However, planning with EVs in medical waste management is more complex than with CEVs given the need to consider recharging.

OPEN IN VIEWER

Table 1.

Recent VRP Studies (2020–2022) on the Waste Collection Problem

Authors	Performance Metrics			FC	AS	Constraints			MO	Methods
	C	E	R			DS	LD	PR
Taslimi et al. ( 21 )	✓	✗	✓	CEV	MW	✗	✓	✗	✗	Exact
Ghannadpour et al. ( 22 )	✓	✓	✓	CEV	MW	✓	✓	✗	✓	EA
Babaee and Aydın ( 23 )	✓	✓	✓	CEV	MW	✓	✗	✗	✓	Scalar
Tirkolaee et al. ( 24 )	✗	✓	✓	CEV	MW	✓	✗	✗	✓	Scalar
Zhao et al. ( 26 )	✓	✗	✓	CEV	MW	✓	✗	✗	✓	Scalar
Govindan et al. ( 25 )	✓	✗	✓	CEV	MW	✓	✗	✗	✓	Scalar
Erdem ( 7 )	✓	✗	✓	EV	MW	✗	✗	✓	✓	Scalar
Erdinç et al. ( 13 )	✗	✓	✗	EV	RW	✗	✓	✗	✗	Exact
Cao et al. ( 14 )	✓	✗	✗	EV	RM	✗	✗	✗	✗	Genetic
Erdem ( 15 )	✓	✗	✗	EV	RM	✗	✗	✓	✗	AVNS
Yang et al. ( 6 )	✓	✓	✗	EV	RM	✗	✗	✗	✗	PS
Masmoudi et al. ( 27 )	✓	✗	✗	EV+ CNG	RM	✓	✓	✓	✗	HTA
Our study	✗	✓	✓	EV	MW	✓	✓	✓	✓	EA+ALNS

Note: AVNS = adaptive variable neighborhood search; VRP = vehicle routing problem; C = cost; E = environmental metric; R = risk; FC = fleet composition; AS = application scenario; DS = trips to disposal site; LD = load-dependent energy consumption; PR = partial recharge; MO = multi-objective; MW = medical waste; RM = residential waste; EV = electric vehicle; CEV = combustion engine vehicle; CNG = vehicle powered by compressed natural gas; EA = evolutionary algorithm; PS = particle swarm; HTA = hybrid threshold acceptance algorithm; Scalar = scalarization-based approaches (e.g., weighted sum method, goal programming).

Electric Vehicle Routing Problem and its Application in Waste Collection

With the growing awareness of the environmental challenge, EVRP has gained more attention in view of EVs’ environmental benefits. Unlike traditional VRP, the EVRP considers driving range, charging strategy, battery capacity, charging duration, and location of recharging stations. Conrad and Figliozzi ( 28 ) were pioneers in introducing EVRP with constant recharge time. Building on this, Macrina et al. ( 29 ) expanded the EVRP to include time windows while using a heterogeneous fleet of CEVs and EVs. More recently, Masmoudi et al. ( 27 ) explored a similar problem using a homogeneous fleet of plug-in hybrid electric refuse vehicles operated by electricity and compressed natural gas (CNG).

Existing research on EVRP can be classified into two categories according to the charging policy: full and partial recharge strategies ( 30 ). The former assumes complete battery replenishment whenever a vehicle reaches a charge station (CS), while the latter allows batteries to recharge up to any level ( 31 ), which is more realistic in real cases. On the issue of energy consumption, some early EVRP studies assumed it was a linear function of travel distance, regardless of the load carried ( 32 ). This assumption is impractical, as consideration of load in discharge calculations can significantly affect routes, particularly when heavy cargo is unloaded early to save energy. Recent research by Vincent et al. ( 33 ), and Masmoudi et al. ( 27 ) considered more realistic energy consumption and partial charging policies. In these studies, energy consumption depends on the weight of the load carried, and the state of charge (SoC) is treated as a decision variable.

Although EVRP has been extensively studied, its application in waste management remains relatively unexplored. Erdinç et al. ( 13 ) first developed a realistic MILP model for the electric garbage truck routing problem using an exact method and real-world data from Istanbul’s Bakirkoy Municipality but they overlooked recharging decisions. Cao et al. ( 14 ) examined an EVRP in a two-echelon recycle network without considering the charging process and proposed a hybrid genetic algorithm combined with large neighborhood search operators. Their findings suggested that using heterogeneous fleets can reduce costs and balance economic, environmental, and resource aspects in recycling systems.

Later, Erdem ( 15 ) introduced an integrated MIP model for WCEVRP with distance-based energy consumption, and developed an algorithm based on adaptive variable neighborhood search (AVNS) in solving this issue. Yang et al. ( 6 ) investigated sorted waste collection with multi-compartment EVs, taking into account the uncertainty of waste generation rates. They introduced a chance-constrained model and proposed a diversity-enhanced particle swarm optimization algorithm to address it. Recently, Masmoudi et al. ( 27 ) examined a single-objective WCVRP using a fleet of vehicles driven by electricity and CNG. They proposed a hybrid threshold acceptance (HTA) algorithm, demonstrating strong performance compared with cutting-edge EVRP algorithms.

Multi-Objective Optimization Methods

Research on multi-objective optimization problems (MOPs) began in the 1970s and, since then, numerous techniques have been developed in solving MO-VRPs ( 34 ). These methods can be divided into two classes: scalar and Pareto methods.

Scalar methods involve mathematical transformations to convert MOPs into single-objective optimization problems. Recently, scalar methods such as the weighted sum approach ( 7 ), goal programming method ( 23 – 25 ), and ϵ -constraint method ( 26 ) have been applied in solving MWVRP. However, these scalar methods have notable drawbacks. For instance, setting weights based on objective importance is challenging, and these methods often produce poorly distributed solutions. Moreover, they can only find Pareto optimal solutions within convex regions ( 35 ), rendering scalar methods potentially ineffective for certain multi-objective optimization problems.

On the other hand, Pareto methods, which use the concept of dominance, have gained popularity in addressing real-world MOPs ( 36 ). This is a result of their superior heuristic optimum-seeking ability and capability to identify a large set of non-dominated points within a reasonable computational time.

In MO-VRPs, Pareto-dominance is often used in conjunction with an evolutionary algorithm (EA) framework, such as NSGA-II ( 37 ), MOEA/D ( 38 ), and SPEA2 ( 39 ). Most EAs in both research and application areas follow a similar framework, using a search operator based on Pareto domination and an iterative reproduction operator ( 40 ). The most used search method is genetic algorithm ( 34 ), other methods such as ant colony ( 41 ), adapted large neighborhood search (ALNS) ( 42 ), GRASP ( 43 ), Tabu search ( 44 ), and particle swarm ( 45 ) have also been combined with Pareto-dominance with remarkable results.

Substantial evidence indicates that incorporating local search mechanisms into EAs significantly improves their performance in combinatorial optimization problems ( 36 ). To our knowledge, no previous studies have used a decomposition-based EA with an ALNS search mechanism to address our specific problem. This research gap motivates our study to develop a novel Pareto-based hybrid algorithm, MOEA/D-ALNS, which combines the MOEA/D framework and ALNS for solving the MO-MWEVRP.

Literature Gap

While there is substantial research related to the waste collection EVRP, there appears to be a limited number of studies that specifically tackle the MO-EVRP within the specific context of medical waste collection and transportation. This subject is particularly crucial given the significant public health and environmental challenges it presents. Our study aims to bridge this gap by concurrently considering routing, waste disposal, and recharging decisions, with special attention to factors such as load-dependent energy consumption, partial recharge policy, and infection risk. The risk function used in this study is determined by variables such as vehicle load, population exposure between nodes, and transmission probability.

Additionally, as shown in Table 1, despite the extensive research on multi-objective WCVRP, which predominantly used scalar methods, a notable research gap exists in the application of evolutionary algorithms to solve the MWEVRP.

Our exploration in this area not only expands the existing MWEVRP literature, but also provides practical solutions to a more complex waste collection routing problem.

Energy Consumption

The energy consumption of an EV is influenced by various factors, including travel distance, mass, vehicle speed, engine efficiency, and so forth. To compute the EV’s energy consumption, we first need to determine the mechanical power requirement, denoted as ME . The computation formula proposed by Demir et al. ( 46 ) can be used to measure this for an EV traveling between node I and j:

M E ij ≈ α ij ( W + L ij ) d ij + β v ij 2 d ij

(1)

where

α ij = a + g sin θ + g C r cos θ is an arc specific coefficient,

β = 0 . 5 C d A ρ is a vehicle coefficient,

W represents the empty vehicle weight (kilograms),

L ij is the load carried by the vehicle on this arc (kilograms), and

d ij denotes the arc ( i , j ) length.

The parameters a , g , C r , C d , A , ρ , and θ stand for the acceleration (assumed to be zero), gravitational constant ( 9 . 81 meter / secon d 2 ) , rolling resistance coefficient (0.01), aerodynamic drag coefficient (0.7), vehicle frontal surface area ( 3 . 912 mete r 2 ) , air density ( 1 . 2041 kilogram / mete r 3 ) , and road gradient (assumed to be zero), respectively. The result is measured in joules ( j = kg m 2 / S 2 ) and can be converted to kilowatt-hours ( kWh ) .

Taking into account the motor efficiency of an EV, denoted as ef f m , the mechanical energy demanded by the EV for each arc is translated into the required electric energy for the vehicle’s motor. This conversion is essential for accurately estimating the amount of battery power needed to supply the corresponding electric energy, considering the battery’s discharge efficiency, represented as eff d . The total amount of battery energy, defined as E ij , required to traverse arc ( i , j ) can be computed as below:

E ij = ef f d · ef f m · M E ij = ef f d · ef f m · [ α ij ( W + L ij ) d ij + β v ij 2 d ij ]

(2)

The values ef f m = 1 . 25 and ef f d = 1 . 11 are proposed by Feng and Figliozzi ( 47 ) and Murakami ( 48 ) in their respective research. For simplicity in the calculation, it is assumed that the travel speed is constant and both the motor efficiency and discharge efficiency are set to 1. Equation 2 can then be converted to the following equation suggested by Rastani and Çatay ( 49 ) for computing energy consumption.

E ij ≈ ( h + w L ij ) d ij

(3)

where h represents the energy consumption per unit distance for an unladen vehicle, while w denotes the extra energy required to transport one unit of cargo.

Mathematical Model

The routing problem in an urban area can be described using a complete and directed network represented as G = ( N , A ) , where N symbolizes the set of nodes (depot, demand location, charge station, and disposal facility) and A stands for the set of arcs connecting these nodes. In this network, i ∈ I represents demand nodes, k ∈ K signifies the vehicle set, and F denotes the set of disposal facility. Furthermore, the set J is comprised of j charging stations.

To allow multiple visits to the same CS, copies of recharging stations are introduced, forming the set J ′ of dummy vertices. The recharging stations j ∈ J ∪ J ′ are defined with a constant charging speed g. Since the disposal facility may be visited more than once, copies of disposal node are generated, resulting in the set F ′ of dummy nodes.

The depot serves as a charge station, and thus, it is included in the set J. By duplicating the depot, the starting and ending nodes are represented by 0 and 0 ′ , respectively. For each arc ( i , j ) ∈ A , a travel time t ij and an associated distance d ij are assigned. This model assumes that charging stations and disposal facilities operate continuously throughout the day.

To better reflect real-world scenarios, our model adopts the partial recharging (PR) policy over the traditional full recharging method ( 6 , 50 ). Figure 1 shows an example with six customers, two CSs, a waste disposal facility, and a depot. The percentages on the routes indicate the vehicle’s battery state of charge on its arrival at a specific location or after recharge at a station. By allowing vehicles to recharge only the necessary energy for their next trip, this policy not only reduces recharge time but also facilitates adherence to the maximum duration constraint. Furthermore, minimizing unnecessary recharging can open up additional customer service time windows thanks to reduced waiting times at charging stations, potentially enhancing the solution quality ( 51 , 52 ). It is noteworthy that under the PR scheme, the recharge quantity, as indicated in Table 2, becomes a continuous decision variable, introducing significant complexity to the problem.

Minimize Obj 1 = ∑ ( i , j ) ∈ A ∑ k ∈ K H ijk x ijk

(4)

Minimize Obj 2 = ∑ ( i , j ) ∈ A ∑ k ∈ K P o p ij L ik R I ij x ijk

(5)

OPEN IN VIEWER

Figure 1.

An illustrative example of a solution for the medical waste electric vehicle routing problem (MWEVRP).

Subject to

H ijk = ( h + w L ik ) d ij ∀ ( i , j ) ∈ A ; i ≠ j , k ∈ K

(6)

∑ j ∈ N j ≠ 0 ′ x jik − ∑ j ∈ N j ≠ 0 x ijk = 0 , ∀ i ∈ N / { 0 , 0 ′ } , k ∈ K

(7)

∑ k ∈ K ∑ j ∈ N j ≠ i x ijk = 1 ∀ i ∈ I

(8)

∑ k ∈ K ∑ j ∈ N j ≠ i x ijk ≤ 1 ∀ i ∈ J ∪ J ′ ∪ F ∪ F ′

(9)

x i 0 ′ k = 0 ∀ i ∈ I ∪ { 0 } , k ∈ K

(10)

∑ j ∈ N j ≠ 0 x 0 jk = 1 ∀ k ∈ K

(11)

∑ i ∈ N i ≠ 0 x i 0 ′ k = 1 ∀ k ∈ K

(12)

L ik = 0 ∀ i ∈ { 0 , 0 ′ } ∪ F ∪ F ′ , k ∈ K

(13)

0 ≤ L ik ≤ C ∀ i ∈ I ∪ J ∪ J ′ , k ∈ K

(14)

L ik + D j ≤ L jk + M ( 1 − x ijk ) ∀ i ∈ N , j ∈ N / F ∪ F ′ ∪ { 0 , 0 ′ } , i ≠ j , k ∈ K

(15)

w ik + ( t ij + s i ) x ijk ≤ w jk + M ( 1 − x ijk ) ∀ i ∈ N / J ∪ J ′ , j ∈ N , i ≠ j , k ∈ K

(16)

w ik + ( t ij + g p ik ) x ijk ≤ w jk + M ( 1 − x ijk ) ∀ i ∈ J ∪ J ′ , j ∈ N / J ∪ J ′ , i ≠ j , k ∈ K

(17)

a i ≤ w ik ≤ b i ∀ i ∈ N , k ∈ K

(18)

0 ≤ w ik ≤ T max ∀ i ∈ N , k ∈ K

(19)

f 0 k = Q ∀ k ∈ K

(20)

0 ≤ f jk ≤ f ik − H ijk + M ( 1 − x ijk ) ∀ i ∈ N / J ∪ J ′ , j ∈ N , i ≠ j , k ∈ K

(21)

0 ≤ f jk ≤ f ik + p ik − H ijk + M ( 1 − x ijk ) ∀ i ∈ J ∪ J ′ , j ∈ N , k ∈ K

(22)

0 ≤ f jk ≤ Q ∀ j ∈ N , k ∈ K

(23)

0 ≤ p ik ≤ Q − f ik + Q ( 1 − x ijk ) ∀ i ∈ J ∪ J ′ , j ∈ N , k ∈ K

(24)

x ijk ∈ { 0 , 1 } , H ijk ≥ 0 , w ik ≥ 0 , p ik ≥ 0 , L ik ≥ 0 , f ik ≥ 0 ∀ i ∈ N , j ∈ N , i ≠ j , k ∈ K

(25)

The objective functions 4 and 5 aim to minimize the overall energy consumption of EVs and the associated infection risk, which is measured by the number of exposed individuals on the road, the vehicle load, and the probability of disease transmission. Equation 6 computes the energy consumption by EV on a given arc. Constraint 7 ensures the input and output flow balance for nodes except the depot. Constraint 8 dictates that each demand node is served once by an EV. Constraint 9 permits vehicles to visit the recharge station and disposal facility multiple times if necessary. Constraints 10 to 12 stipulate that vehicles begin and end their routes at the depot, and the final point visited before returning to the depot is either the disposal facility or the CS. Constraint 13 indicates that the vehicle leaves the depot or the disposal site with an empty load. Constraints 14 and 15 describe the capacity restrictions for vehicles and the load balance between two consecutively visited points, respectively. Constraints 16 and 17 track the arrival time of EVs at a node. Constraint 18 and 19 set the time window for demand nodes and the latest return time for EVs to the depot. Equation 20 indicates that the EV departs the depot with a full battery. Constraints 21 and 22 monitor the SoC for each vehicle. Constraint 23 limits the SoC of an EV to be less than the maximum battery capacity. Constraint 24 accounts for partial battery recharging. Constraint 25 defines the types of variables.

OPEN IN VIEWER

Table 2.

Mathematical Notation

Parameters
d ij	Distance between nodes
D i	Customer demand ( kg )
C	Vehicle capacity ( kg )
P o p ij	Population exposed on arc ( i , j )
g	Recharge speed ( kWh / hour )
Q	Battery capacity ( kWh )
[ a i , b i ]	Time window of demand node i
M	A large number
s i	Service time at node i
R I ij	Disease transmission probability
T max	Maximum working time
Variables
x ijk	Binary variable, which equals 1 when vehicle k travels from node i to j , and 0 otherwise
H ijk	The amount of energy consumed by electric vehicle (EV) k on the arc (i, j)
w ik	The time when vehicle k arrives at point i
p ik	The amount of energy replenished for EV k at recharging station i ∈ J ∪ J ′
L ik	Loading weight when vehicle k departs from node i
f ik	SoC when EV k arrives at node i

General Framework

Algorithm 1 outlines the pseudo-code for the developed MOEA/D-ALNS. The first step defines the algorithm’s inputs, including: i) S : the total number of subproblems produced by the process of decomposition; ii) N : the number of neighbors for each subproblem; iii) η max and η noi max : the maximum number of iterations and non-improvement iterations; and iv) ALNS-related parameters ( Ω − , Ω + , p − , p + , q , α and T ) are defined.

OPEN IN VIEWER

Algorithm 1 MOEA / D − ALNS

Input: Subproblem ( S ) , neighbor num ( N ) , maximum iteration num ( η max ) , maximum number of non-improvement iterations ( η noi max ) , destroy and repair operators ( Ω − , Ω + ) , selection probability ( p − and p + ) , num of maximum removal nodes ( q ) , Temperature ( T ) , cooling rate ( α ) Output: Set of non-dominated solutions ( D ) 01: W ← Initialize weight ( S ) #Scheffé’s method 02: G ← Neighbors generation ( W , N , S ) 03: D ← Ø #Record Non dominant solution 04: Obj_value ← Ø #Record the objective value of solution 05: T ← T 0 06: for i in S do 07:  Po p i ← Generate initial solution ( ) 08:  so l i ← Apply the ISD ( so l i ) and IRS ( so l i ) 09: end 10 z * ← update ( so l i ) 11: while n < η max or n < η noi max do 12: for i in S do 13:   Choose operators d ∈ Ω − and r ∈ Ω + 14:   so l i ← r ( d ( so l i , Ω − , Ω + , q , p − , p + ) ) 15:   so l i ← Local search ( so l i ) 16:   so l i , Po p i ← ISD ( so l i ) and IRS ( so l i ) 17:   z * ← update ( so l i , z * ) 18:   Update scores and weights for Ω − and Ω + 19:  for neighbor j in G i do 20:   if g te ( so l i | w j , z * ) ≤ g te ( so l j | w j , z * ) 21:      Po p j ← Po p i 22:    end 23:  end 24:   Obj_ value ← Obj_ value ∪ Obj_ value { so l i } 25: end 26:  Update selection probability p − , p + 27:  T ← T * α 28: end 29: for i in S do 30:  D ← Non_dominant_sort ( D , Obj _ value ) 31: end 32: return D

Next, the MWEVRP is decomposed into S scalar subproblems using the weight set W and Scheffé’s method, adhering to the condition: | W | = W 1 + W 2 + ⋯ + W r = 1 , where r denotes the number of objectives. This approach is also used in the classic MOEA/D ( 59 ). The set G i stores the N nearest neighbors for each subproblem i .

To generate an initial individual solution for each population ( Po p i ) while maintaining Pop’s diversity, customers are randomly assigned to different EVs until all customers have been allocated. Subsequently, two improvement operators (IDS and IRS) for waste disposal and recharge constraints are employed to ensure solution feasibility. The generated solution is assessed for uniqueness, and if unique, it is added to the population. This approach creates an initial population of size S. After generating the population, the best solution for each objective forms a vector z * .

During each iteration, the algorithm performs the following steps:

On meeting the termination condition of the loop ( η max and η noi max ) , the algorithm returns the set of Pareto solutions based on the non-dominated sorting principle.

min v ∈ Ω g t e ( v | λ , z * ) = max 1 ≤ i ≤ n { λ i | f i ( v ) − z i * | }

(27)

Tchebycheff Approach Aggregation Function

The weighted Tchebycheff approach can handle problems with non-convex Pareto fronts ( 59 ). A scalar optimization subproblem is formulated as Equation 27.

The ideal reference point is defined as z * = ( z 1 * , z 2 * , … , z n * ) (if z i * < min { f i ( v ) } | v ∈ γ , i = 1 , … , n ) , and the weight vector λ = ( λ 1 , λ 2 , … , λ n ) , so that ∑ i = 1 n λ i = 1 , λ i ≥ 0 is true. The solution space is defined by the variable v , and n stands for the number of objective functions.

Adaptive Large Neighborhood Search

ALNS has demonstrated remarkable performance in addressing diverse VRPs ( 33 ). To generate the neighborhood of a given solution, six destroy and four repair operations are performed. Following each iteration, a scoring mechanism evaluates and updates the past performance of operators. Operators are chosen according to their scores, and high-scoring operators are favored for subsequent iterations. A newly generated solution is accepted based on a predefined acceptance criterion. Unlike the single-objective ALNS mechanism, this study employs the TCH aggregation function to evaluate solution quality.

Adaptive Mechanism

The adaptive mechanism modifies the probability P of selecting each destroy and repair operator based on its past performance. During each iteration of the ALNS, an operator’s score increases by σ a for a new improving solution, σ b for a new deteriorating but accepted solution according to the simulated annealing (SA) criterion, and σ c for no improvement. Initially, probabilities for six destroy and four repair operators are set. During the algorithmic procedure, the weight w j are updated as:

w j n + 1 = ( 1 − r ) w j n + r π j n

(28)

This Equation 28 uses the parameter r to balance the importance of the current performance and the past performance of the operator, where r ranges from 0 to 1, and π j n is defined as:

π j n = score of operator j after n th iteration total times operator j selected after n th iteration

(29)

Specifically, the probabilities of destroy operators ( p − ) and repair operators ( p + ) are periodically adjusted, the corresponding probabilities to choose operators are derived by: w j ∑ j w j . Once the weight is updated, all variables in Equation 29 are reset to 0.

Simulated Annealing Acceptance

To avoid local optimality, the SA is employed. In each iteration, a probability is computed using Equation 30 to accept a new, worse solution. The new solution is accepted if the probability exceeds μ , where μ is a random number ranging from 0 to 1. The temperature T is initially set and decreased after each iteration by multiplying it with a cooling parameter α ( 0 < α < 1 ) .

probability = e − ( g te ( so l new | w j , z * ) − g te ( so l current | w j , z * ) ) T

(30)

Method for Updating Non-Dominated Solutions

This research integrates multiple objectives into the optimization process to generate solutions. However, since these objectives may conflict, finding an optimal solution that minimizes all of them concurrently is not possible. Therefore, we focus on identifying the Pareto frontier, which comprises solutions that are not dominated by any other solution. To update the set of non-dominated solutions (D), the non-dominant sort method presented in Cota et al. ( 54 ) is employed. A solution vector A = [ x 1 , x 2 , … , x n ] dominates another vector B = [ y 1 , y 2 , … , y n ] if and only if:

∀ i ∈ [ 1 , … n ] : A i ≤ B i and ∃ i ∈ [ 1 , … , n ] : A i < B i

(31)

Parameter Tuning

Parameter calibration is achieved through empirical tests on five instances (Table 3). Each factor is evaluated individually while keeping other parameters constant. We set the stopping criterion for the MOEA/D, MOEA/D-LNS and MOEA/D-ALNS algorithms at 20 generations ( η max = 20 ) , as empirical tests within a range of 10 to 50 generations showed that this criterion was essential for convergence. The parameter S denotes the number of subproblems created in all algorithms, which influence the granularity of weights during the search process. Empirical tests for S values between 5 and 100 revealed that setting S to 20 for small instance, 30 for medium instance, and 50 for large-scale problems resulted in fast solution convergence. Finally, the parameter N is tested between 2 and 20, and optimal values of 6, 8, and 10 were identified for small-, medium-, and large-scale instance, respectively, to achieve good results.

OPEN IN VIEWER

Table 3.

Data About the Test Instances

Instance	D	CSs	V	Pop	RI
c202C10	10	5	2	(50,100)	(0.001,0.002)
rc202C15	15	5	2	(50,100)	(0.001,0.002)
c201_ 21_25	25	21	3	(50,100)	(0.001,0.002)
rc201_21_50	50	21	5	(50,100)	(0.001,0.002)
rc201_21_100	100	21	8	(50,100)	(0.001,0.002)

Note: D = no. of demand nodes; CS = no. of charge stations; V = no. of vehicles; Pop = population exposure between nodes; RI = infection probability between nodes.

In the NSGA-II algorithm, both the crossover and mutation rates are two pivotal parameters that significantly influence the algorithm’s performance ( 63 ). Testing these rates is crucial to determine if a specific combination yields satisfactory results. For this purpose, various values of the crossover and mutation rates were used to obtain solutions for an instance with 15 customers. Following this, BEC, BR, HV, and CPU values are calculated to evaluate and compare the quality of the solutions. This evaluation aids in identifying which combination of parameter values produces the most optimal solution.

Table 4 displays the metrics for each tested parameter combination. Based on these results, the optimal values determined for the crossover and mutation rates are 0.9 and 0.1, respectively. This combination yielded the best performance for BEC and BR, and also demonstrated favorable results in HV. To ensure a fair comparison, the parameters of population size ( S ) and maximum iteration ( η max ) used in NSGA-II are set identically to three other algorithms (MOEA/D, MOEA/D-LNS, and MOEA/D-ALNS).

OPEN IN VIEWER

Table 4.

Parameter Test for NSGA-II.

Crossover	Mutation	BEC	BR	HV	CPU
0.9	0.3	488.2	129.5	11.6	103.4
	0.2	489	132	13	96
	0.1	464	129	11.5	93.3
0.8	0.3	493	130	11.8	71.6
	0.2	491.5	137.5	9.81	76.1
	0.1	514.2	132.5	9.83	73.5
0.5	0.3	525.7	135.3	9	50.6
	0.2	545.8	139.4	11.4	48 . 2
	0.1	536.6	136.6	9	48 . 2
0.2	0.3	inf	inf	na	na
	0.2	inf	inf	na	na
	0.1	inf	inf	na	na

Note: NGSA = non-dominated sorting genetic algorithm; BEC = best energy consumption; BR = best risk; CPU = computational time; HV = hypervolume; Inf = infeasible solution; na = not applicable.

Results and Discussion

In this section, the performance of MOEA/D-ALNS and the enhancement strategies for MOEA/D are evaluated. Specifically, MOEA/D-ALNS is compared with the following algorithms: i) the NSGA-II enhanced with specialized crossover and mutation operators, and ii) the modified MOEA/D and MOEA/D-ALNS without the adaptive selection strategy.

Comparison between MOEA/D-ALNS and NSGA-II

In this subsection, MOEA/D-ALNS is compared with the NSGA-II algorithm, with a focus on the objectives BEC and BR. Results for small and medium instances are shown in Table 5. For the large-scale instance, when both algorithms used identical values for population size and maximum iteration, NSGA-II failed to generate viable solutions. Consequently, these results were omitted from the table. Table 5 reveals that MOEA/D-ALNS consistently surpasses the NSGA-II algorithm for both objectives, particularly for medium-sized instances. For example, MOEA/D-ALNS achieves an energy saving of approximately 37 % for the objective of minimizing total energy consumption. Additionally, it reduces the total risk by 9 % compared with the NSGA-II algorithm in instances with 50 customers. With regard to computation time, NSGA-II performs better. However, the time difference between these two algorithms narrows with increasing problem size.

OPEN IN VIEWER

Table 5.

Comparisons between MOEA/D-ALNS and NSGA-II on Objectives

Instance	Minimizing the energy consumption		Minimizing the risk		Computational time
	NSGA-II	MOEA/D-ALNS	NSGA-II	MOEA/D-ALNS	NSGA-II	MOEA/D-ALNS
c202C10	254	254	79	79	15 . 52	18.43
rc202C15	484	477	132.14	129 . 29	102 . 21	140.03
c201_21_25	696.8	548 . 1	227.1	220 . 2	172 . 03	202.08
rc201_21_50	1,775.4	1112 . 2	507.8	461 . 4	641 . 83	643.3

Note: MOEA/D = multi-objective evolutionary algorithm using decomposition; ALNS = adaptive large neighborhood search; NSGA = non-dominated sorting genetic algorithm.

Comparison of MOEA/D, MOEA/D-LNS, and MOEA/D-ALNS Results

A comparison of the results in Tables 6 and 7 reveals that for energy consumption, the three algorithms show no significant differences in small data sets with 10 and 15 demand nodes. However, as the problem dimension increases, MOEA/D-ALNS consistently outperforms the others, being the most energy-efficient with the lowest average energy consumption (921.61), followed by MOEA/D-LNS (976.01) and MOEA/D (1,097.08).

OPEN IN VIEWER

Table 6.

Comparisons of MOEA/D-Based Algorithms Based on BEC, BR, and CPU

Instance	MOEA/D			MOEA/D-LNS			MOEA/D-ALNS
	BEC	BR	CPU	BEC	BR	CPU	BEC	BR	CPU
c202C10	254	82	6.13	254	82	9.75	254	82	9.78
rc202C15	524.2	169.7	54.73	535.1	171.1	52.07	543.7	174.3	68.29
c201_21_25	507.29	220.43	116.13	541.66	223.36	182.35	549.83	219.64	149.30
rc201_21_50	1,258.7	494	441.39	1,165.2	501.7	478.25	1,065.6	496.6	596.12
rc201_21_100	2,941.2	1,306.2	3,376.51	2,384.1	1,275.7	3,780.75	2,194.9	1,263.2	3,839.45
Average	1,097.08	454.47	798 . 98	976.21	450.77	900.63	921 . 61	447 . 15	932.59

Note: BEC = best energy consumption; BR = best risk; CPU = computational time; NPS = number of Pareto solutions; HV = hypervolume; MOEA = multi-objective evolutionary algorithm; MOEA/D = multi-objective evolutionary algorithm using decomposition; LNS = large neighborhood search; ALNS = adaptive large neighborhood search.

OPEN IN VIEWER

Table 7.

Comparisons of MOEA/D-Based Algorithms Based on NPS, HV, and SP

Instance	MOEA/D			MOEA/D-LNS			MOEA/D-ALNS
	NPS	HV	SP	NPS	HV	SP	NPS	HV	SP
c202C10	2	0	0	2	0	0	2	0	0
rc202C15	2.6	8.7	0.2	2.1	8.1	0.34	2.8	6.2	0.20
c201_21_25	3.71	9.43	0.23	3	10.14	0.22	2.79	9.14	0.17
rc201_21_50	6.1	11.7	0.13	4.8	12.3	0.23	5.1	13.5	0.23
rc201_21_100	5.6	12.3	0.18	5.0	15.4	0.25	5.3	15.9	0.20
Average	4 . 00	10.53	0 . 19	3.38	11 . 48	0.26	3.60	11.19	0 . 19

Note: NPS = number of Pareto solutions; HV = hypervolume; SP = spacing; MOEA = multi-objective evolutionary algorithm; MOEA/D = multi-objective evolutionary algorithm using decomposition; LNS = large neighborhood search; ALNS = adaptive large neighborhood search.

The average values for HV and SP were calculated without using data from the c202C10 instance.

With respect to risk metric, the three algorithms perform similarly, with MOEA/D-ALNS having a minor advantage in the large-scale instance. MOEA/D demonstrates superior NPS and CPU performance, indicating its effectiveness in generating diverse Pareto solutions and its suitability for situations when fast computation is a priority. In the HV comparison, MOEA/D-LNS outperforms the other two algorithms, while MOEA/D and MOEA/D-ALNS have similar performance in SP, which is slightly better than MOEA/D-LNS.

The radar plot in Figure 3 displays the normalized means of the metrics for these three algorithms, using a distance-based normalization approach ( 64 ). The greatest performance in this scenario obtains one point, while the worst performance obtains 0. The computational results reveal that MOEA/D-ALNS is ideal for obtaining the most energy- and risk-efficient solutions, while MOEA/D is tailored for prompt solution generation.

OPEN IN VIEWER

Figure 3.

The radar plot for three algorithms.

The three algorithms are further evaluated using the simple additive method (SAM) with six equally weighted performance metrics. These metrics are normalized using a distance-based technique ( 64 ). The ranking of these three EAs can be found in Table 8, with the MOEA/D-ALNS algorithm receiving the highest score of 0.68. The result indicates that the MOEA/D-ALNS algorithm outperformed and dominated the other two EAs. The MOEA/D and MOEA/D-LNS algorithms rank second and third, respectively.

OPEN IN VIEWER

Table 8.

The Ranking of Algorithms Based on SAM

Algorithm	Indicator						Score
	BEC	BR	NPS	HV	SP	CPU
MOEA/D	1,097.08	454.47	4.00	10.53	0.19	798.98	0.50
MOEA/D-LNS	976.01	450.57	3.38	11.48	0.26	900.46	0.41
MOEA/D-ALNS	921.61	447.15	3.60	11.19	0.19	932.58	0 . 68

Note: SAM = simple additive method; BEC = best energy consumption; BR = best risk; CPU = computational time; NPS = number of Pareto solutions; HV = hypervolume; SP = spacing; MOEA = multi-objective evolutionary algorithm; MOEA/D = multi-objective evolutionary algorithm using decomposition; LNS = large neighborhood search; ALNS = adaptive large neighborhood search.

Notably, our results closely align with a recent experiment conducted by Wang et al. ( 63 ) which focused on the multi-objective perishable product VRP. Unlike our approach, their research developed an algorithm that combined the NSGA-II framework with ALNS operators. This new approach, termed NSGA-II-ALNS, was compared with the basic NSGA-II algorithm and the NSGA-II-LNS algorithm. Their findings demonstrated superior overall performance of NSGA-II-ALNS over both the basic NSGA-II and the NSGA-II-LNS. When comparing the performance of the LNS algorithm under the NSGA-II and MOEA/D frameworks with the basic NSGA-II and MOEA/D, respectively, both our study and that conducted by Wang et al. ( 63 ) presented consistent results. In particular, the inclusion of the LNS algorithm led to an improvement in the objective function and an increase in HV. This consistency across different studies further validates the notion that integrating ALNS operators into the evolutionary framework, along with other tailored improvement strategies, can significantly enhance the quality of derived solutions.

Figure 4 shows the Pareto fronts for all algorithms, illustrating the dispersion of Pareto points in the solution area. An inverse correlation between the two objective functions is observed. This implies that reducing energy consumption can lead to a rise in the risk function, and vice versa. In practical terms, the optimal solution concerning the risk objective function might come at a high energy consumption. Conversely, the solution with the lowest energy consumption is likely to involve a high-risk route.

OPEN IN VIEWER

Figure 4.

The evolutionary-based algorithms’ Pareto fronts: (a) instance: rc201_21_50; (b) instance: rc_21_100.

Sensitivity Analysis and Implications

This section investigates the effects of time windows on the objective function by removing them. Additionally, a sensitivity analysis is performed to evaluate how changes in battery capacity affect the objective functions. If the impact is minor, a non-fully charged policy can be implemented to preserve battery health. Since the previous tests revealed that the MOEA/D-ALNS algorithm achieved better solution quality than other algorithms, it is used exclusively for all tests in this section.

The first experiment examines the influence of time windows on the solution. Table 9 demonstrates that when comparing scenarios with and without time windows, the average energy consumption decreased by 35.56%, 26.87%, 40.86%, and 51.21%, while the average risk values reduced by 29.43%, 17.32%, 12.42%, and 15.95% for the instances with 15 , 25 , 50 , and 100 demand nodes, respectively. Furthermore, CPU significantly decreased in all experimental scenarios, with reductions of 12.44%, 37.08%, 56.47%, and 58.64% for these four instances, respectively.

OPEN IN VIEWER

Table 9.

Sensitivity Analysis Results for Time Windows

Instance	BEC	BEC without TW	% Decrease	BR	BR without TW	% Decrease	CPU	CPU without TW	% Decrease
rc202C15	543.7	350.4	−35.56 %	174.3	123	−29.43 %	68.29	59.79	−12.44 %
c201_21_25	549.83	402.1	−26.87 %	219.64	181.6	−17.32 %	149.30	94.03	−37.08 %
rc201_21_50	1,065.6	630	−40.86 %	496.6	435	−12.42 %	596.12	259.23	−56.47 %
rc201_21_100	2,194.9	1,070.5	−51.21 %	1,263.2	1,061.7	−15.95 %	3,839.4	1,588.71	−58.64 %

Note: BEC = best energy consumption; BR = best risk; CPU = computational time; TW = time window.

Relaxing time-window constraints in the MWEVRP yields substantial improvements in energy consumption, risk value, and computational efficiency. Particularly, when the number of demand nodes increases to 50 or more, the CPU required to solve the problem reduces by more than half. As depicted in Figure 5, these advances are likely a result of the increased flexibility in visiting nodes and better usage of charging stations. Without time-window constraints, vehicles can minimize idle time and continuously move between nodes, thus establishing more energy-efficient and safer routes. Additionally, the removal of customer service time constraints allows the algorithm to explore the solution space more efficiently.

OPEN IN VIEWER

Figure 5.

Time windows impact on routes (c201_21_25): (a) included; (b) excluded.

The imposition of time-window constraints in medical waste management primarily arises from pressing public health concerns ( 12 ). These constraints ensure that hazardous wastes are disposed of within strict timeframes, thus mitigating associated health threats. Normally, medical waste collection is synchronized with the operational hours of health care institutions. However, unforeseen events, such as disease outbreaks, can trigger surges in waste volumes. Such unexpected surges challenge the established time constraints, leading to potential backlogs and heightened safety risks.

Given the challenges posed by strict time-window constraints, emerging practices present promising solutions. Notably, not all medical waste carries the same level of risk. Wastes that have undergone initial treatment, or those without harmful bio-agents, can be managed with more flexibility, reducing the need for rigid time constraints ( 65 ). Furthermore, with a rise in health care facilities operating around the clock, off-peak collection through outsourcing is becoming an increasingly viable solution, addressing the logistical challenges of tight time windows ( 23 ). The integration of these practices (such as flexible time management of low-risk waste and off-peak collection) into our research offers a practical path to alleviate time-window constraints, thereby enhancing the efficiency and safety of medical waste management.

The second experiment involves a detailed sensitivity analysis to explore the relationship between the objective functions and potential changes in battery energy levels. Reserved energy levels (REL) of 10 %, 7.5 %, 5 %, and 2.5 % are considered, which indicate that 90 %, 92.5 %, 95 %, and 97.5 % of the energy stored in the battery can be consumed by the EV. The MOEA/D-ALNS algorithm is used to solve the modified problem instances, with the results recorded in Table 10. Table 11 further illustrates the relative percentage deviation, minimum, maximum, and average change rates of objective functions as the REL varies from 10 % to 2.5 %. These indices are designed to identify potential fluctuations and overall trends in energy consumption and risk. The relative percentage deviation is formulated as follows:

R P D i = g i ( R E L 1 ) − g i ( R E L 0 ) g i ( R E L 0 ) × 100 %

(32)

where

RP D i represents the relative percentage deviation for the i th objective function,

g i denotes the value of the i th objective function at a specific REL, and

g i ( RE L 0 ) and g i ( RE L 1 ) represent the values of the i th objective function at the initial and final REL, respectively.

OPEN IN VIEWER

Table 10.

Average Results of Problem Instances with Various Reserved Energy Level

Instance	Obj	Reserved energy level
		10 %	7 . 5 %	5 %	2 . 5 %	0 %
rc202C15	BEC	1,383	506 . 16	534.67	510.29	543.7
	BR	1,027.57	191.33	193.5	191.29	174 . 3
c201_21_25	BEC	555.625	536 . 83	550.83	588.83	549.83
	BR	234.375	231.67	229.00	231.5	219 . 64
rc201_21_50	BEC	1,173.5	1,117.33	1,112.33	1,090.28	1065 . 6
	BR	507.67	502.67	502.67	495 . 57	496.6

Note: BEC = best energy consumption; BR = best risk.

OPEN IN VIEWER

Table 11.

Relative Percentage Deviation of Objective Functions Between Various Reserved Energy Levels

Instance	Obj.	10 % – 7 . 5 %	7 . 5 % – 5 %	5 % – 2 . 5 %	2 . 5 % – 0 %	Min.	Max.	Avg.
rc202C15	BEC	−63.41%	5.63%	−4.55%	6.54%	4.55%	63.41%	−13.95%
	BR	−81.37%	1.13%	−1.14%	−8.87%	1.13%	81.37%	−22.56%
c201_21_25	BEC	−3.38%	2.61%	6.90%	−6.62%	2.61%	6.90%	−0.12%
	BR	−1.15%	−1.15%	1.09%	−5.13%	1.09%	5.13%	−1.59%
rc201_21_50	BEC	−4.79%	−0.45%	−1.98%	−2.26%	0.45%	4.79%	−2.37%
	BR	−0.98%	0.00%	−1.41%	0.21%	0%	1.41%	−0.55%

Note: BEC = best energy consumption; BR = best risk; Min. = minimum change rate; Max. = maximum change rate; Avg. = average change rate; Obj. = objective.

Our analysis reveals an inconsistent relationship between the objective functions and the REL, despite a noticeable trend of declining objective function values as the REL decreases. This inconsistency seems to be driven by both the size of the considered instance and the inherent characteristics of the proposed solution algorithm.

For instance size, the average change rates, as outlined in Table 11, show that in small instances such as c202C15, where only a few vehicles are deployed, both energy consumption and risk are highly sensitive to shifts in the REL. An increase in REL might necessitate more frequent charging for EVs, thereby leading to a significant surge in both energy consumption and risk. For instance, in the rc202C15 scenario, substantial reductions in energy consumption and risk are observed when the REL shifts from 7.5% to 10%. In contrast, in larger instances like c201_21_25 and rc201_21_50, adjustments in the REL exert only a marginal impact on energy consumption and risk. This may be attributed to the larger fleet size, offering more total battery capacity and flexibility in route planning.

Apart from the instance size, the observed inconsistency may also stem from the inherent characteristics of the solution algorithm being used. This algorithm produces solutions that involve multiple conflicting objectives, aiming to identify an optimal compromise between them. As a result, it can lead to outcomes that are not always consistent. For instance, it could generate a solution with low risk but high energy consumption, or vice versa, resulting in fluctuations in these two objective functions.

Metrics such as the minimum and maximum change rates in Table 11 help illustrate the degree of variability in the solution across different data instances, indicating that as the instance size increases, the variation in the objective function values tends to diminish.

In practice, it is advisable to keep the battery energy level slightly below its full capacity to optimize energy consumption and risk management, while also preserving battery health. Our results suggest that an energy level of 92.5 % seems optimal for small instances with 15 or 25 demand nodes, while a level of 97.5 % might be more suitable for larger instances with 50 demand nodes. Thus, determining the ideal energy level requires a comprehensive evaluation of specific problem instances to attain an optimal balance between energy consumption, risk, and battery health.