AN ATTRACTORS-BASED PARTICLE SWARM OPTIMIZATION FOR MULTIOBJECTIVE CAPACITATED VEHICLE ROUTING PROBLEM

In this paper, a new multiobjective discrete particle swarm algorithm is presented for the Capacitated vehicle routing problem. The binary algorithm integrates particle displacement based on local attractors, a crowding distance as elitism policy and genetic operators. The proposed approach is first implemented on a set of well-known benchmarks for single-objective capacitated vehicle routing problems and compared to results underlined in the literature. The obtained results demonstrate its ability to achieve the main optimization solution and sometimes prove its efficiency from other given techniques. Then, the approach is applied to an academical example for the multiobjective Urban Bus Routing Problem with Route Balancing. Mathematics Subject Classification. 90C29, 90B20. Received February 8, 2021. Accepted August 4, 2021.


Introduction
The vehicle routing problem (VRP) is a combinatorial optimization problem seeking to service a number of customers with a fleet of vehicles. Introduced first by Dantzig and Ramser [8], the VRP is an important problem in the fields of transportation, distribution, and logistics. The problem is concerned with delivering goods to a set of customers with known demands through vehicle routes by minimizing the total route cost. Capacitated vehicle routing problem (CVRP) is the most elementary version of the VRP and seeks a number of routes for number of vehicles with units homogeneous/inhomogeneous capacity to minimize total transportation cost of routes while satisfying the delivery demands of number of customer nodes. Each route must start and finish at the depot and each customer has to be visited exactly once by one vehicle.
Various real cases versions of VRP can be found in literature and details with formulation, and solution methods can be found in [35]. The VRP is formulated as a combination of the bin packing problem (BPP) and the traveling salesman problem (TSP). Since both BPP and TSP are proven to be NP-hard problems [38]; the VRP and its CVRP variant are too. Among those CVRP variants are Vehicle Routing Problem with Simultaneous Pick-up and Delivery (VRP-SPD) [33] and Vehicle Routing Problem with time windows [9,37].
Keywords. Particle swarm optimization, multiobjective optimization, vehicle routing problems, urban bus routing problems, simultaneous pick-up and delivery.
In the VRBP-SPD version, customers require simultaneous collection of goods from their location in addition to distribution of goods to their location and a fleet of vehicles originated in a depot serves customers with pickup and deliveries from/to their locations. Increasing importance of reverse logistics activities make it necessary to determine efficient and effective vehicle routes for simultaneous pick-up and delivery activities. Moreover, when taking into account the costumers waiting time, the vehicle running time and seeking to minimize that time for given reasons (using delay for the delivering goods, preserving some physical and/or medical qualities, reducing the vehicle motor running and the availability time windows of customers), we formulate a multiobjective optimization problem. That is, the new goal is not only to minimize the total transportation cost of route but also minimizing each driver's working time. This multiobjective version of the VRP is named the Vehicle Routing Problem with Route Balancing (VRP-RB) [17]. Numerous real-life applications of transportation problems can be derived from this version. We can cite for instance the Urban School Bus Routing [13,29] which is a challenging problem in urban densely populated areas.
Many natural phenomenon based techniques, for example ant colony optimization (ACOs) and genetic algorithms (GAs), have been developed in solving multi-objective optimization problems such as NSGA-II [10], PAES [20] and SPEA2 [40] in term of genetic algorithms and for Multiobjective Ant colony Optimization, one can refer to [2,26] among others. These population-based algorithms have the ability to explore the different parts of the Pareto front simultaneously [16].
During these last decades, many PSO [18] approaches have been proposed and most of them use the concept of Pareto dominance to classify the solutions [22,27]. Comparisons of different multiobjective algorithms that use such techniques have been outlined in [21,27] among others. These outlined works show that, in general, MOPSO based techniques are able to generate a best set of nondominated solutions close to the true Pareto front compared to GAs based algorithms. Note that in terms of diversity of the nondominated solutions, GA based algorithms can produce best results but are not able to cover the entire Pareto front in all test functions.
However, numerical approximation challenges reside on proposed optimization problems algorithms, such as how intensifying the final computed Pareto front close to the true Pareto known front and how improve diversification of these optimal solution over the whole Pareto front. First, it should be noted that the use of a new displacement of particles [15], which basically is based on local attractors, takes into account more the process of intensification and ensure optimal solution converge to a slongtable point of the true Pareto front. In other hand, the use of crowding distance for the nondominated sorting in [27], the local choice of global best guide of a particle and the personal best guide mutation coupled to a crossover operator may improve considerably the process of diversification.
In other hand, several combinatorial Particle Swarm Optimization (PSO) techniques have been proposed to solve different versions of the Vehicle Routing Problem [1,24,25] and results have shown that PSO based algorithms are able to solve efficiently this kind of combinatorial optimization problem. So, the proposed algorithm extends PSO in solving combinatorial multiobjective capacitated vehicle routing problem with route balancing by changing the classical displacement to the one based on local attractors and incorporating the mechanism of crowding distance computation in the insertion method and the deletion method of the external archive of nondominated solutions. The proposed approach is named the Discrete Multiobjective Particle Swarm Optimization based on Attractors (DMOPSO-A). The attractors based displacement, the crowding distance mechanism together with a mutation based on personal best guides and crossover operators maintain the diversity and the intensity of nondominated solutions in the external archive.
In order to evaluate the performance of the proposed DMOPSO-A approach, several test problems from the literature benchmarks are considered. Results of some known solvers are compared with the solution provided by the proposed approach. And then, an academic version of the urban bus routing problem is solved by the actual approach.
This paper is organized as follows; in Section 2, the mathematical formulation and details of CVRP with route balancing is presented. In Section 3, details of the proposed DMOPSO-A approach are introduced and explained. Numerical experiments and performance validations of the proposed approach are evaluated in Section 4. Finally, Section 5 concludes the paper with future perspectives.

The multiobjective capacitated vehicle routing problem
In this section, the CVRP is formulated as a Mixed-integer linear programming (MILP) model based on model. Let = ( , ) be a graph with = ( 0 , 1 , . . . , ) a set of vectors where 0 is the depot and 1 , . . . , are the customers to be served; and = = {( , ) : ( , ) ̸ = ∈ } a set of arcs, representing a route from customer to customer . The graph is considered to be complete, . ., for each couple ( , ) ̸ = of customers, there exists an arc (a direct route) from customer to customer . Matrices and with entries and represent the length of the arc and the average time for the road cross respectively. Each customer has a given positive delivery , a pick-up quantity and a service duration . At the depot, vehicles are located and each vehicle has a maximal capacity and a maximal service duration . Each vehicle has also a fixed cost and a variable cost per distance unit.
The multiobjective CVRP consists on choosing at most circuits so that (1) Each circuit starts and ends at the depot; (2) Each customer is visited once with exactly one vehicle; (3) The total discharge in a vehicle along a given arc must not exceed the maximal capacity of the vehicle; (4) The total service duration of each vehicle (including the service duration) must not exceed the limit duration of the vehicle; (5) The total cost and the route balancing of the whole routing must be as minimal as possible.
Let define the CVRP decision variables as a binary variable which indicates whether the arc ( , ) is traversed by the vehicle or not. = 1 is the vehicle follows the route ( , ) and = 0 if not; and a variable which represents the weight of vehicle when traversing the route ( , ). The objective functions of this problem are expressed as: and the optimization problem is stated by 3. The binary particle swarm optimization based on attractors 3.1. The discrete particle swarm optimization The initial binary particle swarm optimization algorithm was proposed by Kennedy and Eberhart [19] in 1997 for discrete optimization problems and later other versions are studied in [5,6,28] among others. In these methods, each particle must take "true" or "false" values, "yes" or "no" values, or solution alike which indicate a potential solution. A fitness solution is always introduced in order to evaluate the appropriateness of solutions. Each particle is located in some dimensional space and each element of a particle position can take the binary value of 0 or 1 which can change from 1 to 0 and vise versa. Each particle has also a velocity vector with elements defined in range [− max , max ]. The DPSO velocities are defined in terms of probabilities that a bit will be in one state or the zero state. At the beginning of the algorithm, a number of particles and their velocities are initialized randomly. Then in each iteration step, the algorithm searches on how reaching the optimal solutions based on its predefined fitness function. The global best guide and the personal best guides used to update velocity vectors are chosen in the same manner as for the continuous optimization. Eventually, the position of the particles is updated using velocity vectors.
With this definition, the velocity vector must be restricted in the interval [0, 1]. So, a normalization and control function is introduced to force the real values of velocity vectors to be in [0, 1]. In this study, the sigmoid function [36] is defined by The velocity is updated in the same way as in continuous optimization and is expressed by where is the inertia factor that influences the local and global abilities of the algorithm, , is the velocity of the particle at iteration step , 1 and 2 are weights affecting the cognitive and social factors, respectively.
1 and 2 ∼ (0, 1); , stands for the best value found by particle (pBest) and , denotes the global best found by the entire swarm (gBest). And the new position of the particle is obtained using the equation where is a uniform random number in the range [0, 1].

Displacement based on attractors
Attractors, in PSO, describe the path marked by particle positions which any particle passes by before landing on its final target.  .3), the best personal performance , and the best global performance , of the swarm. So far the particle ( , +1 ; , ) travels by the local positions ( , and , ) before it reaches the next particle position ( , +1 ; , +1 ). The operators which recently inserted in [15] perform the continuous and/or combinatorial optimization. The increment process grants the application a benefit of the convergency and the diversity. The equations of attractors are stated as follows: The constants 1 , 2 , and 3 are arbitrary numbers and 1 , 2 and 3 ∼ (0, 1). However, it is suitable to secure the convergency of the approach. Therefore, the arbitrary constants will be constraint as follows (one set at a time) [15]: Introducing the sigmoid function, the attractors aforementioned are reshaped to correspond as follows. The velocity remains as it is typed above (3.4), however, the intermediate position are formulated so that the following particle position is written where ( , , ) ∼ (0, 1) are random numbers and | | is the absolute value of .

Elitist policy
The diversity is a quality performance which quantify an approach for its abilities to spread particles and have a comparable Pareto front. Hypercubes techniques [7] and Crowding Distance (CD) [23,27] are well known in operations. They aim to select the best particles in order to respect an adequate distribution of the Pareto. The process consists of mining particles nondominated solutions. The following work, The CD gathers to perform this tedious work and it delivers value of density of solutions concerning its surroundings. The obtained values are used to for comparison amongst nondominated solutions, therefore to maintain a well-distributed Pareto front. To illustrate, if a given archive set is full, one of the techniques (hypercubes techniques or CD) is applied to select the best solution between a nondominated solution in the archive set and a new nondominated solution and stored in the archive.
For instance, the nondominanted solutions rest in the archive set; they are arranged in order, from lower to the upper of objective function values. These values are computed for all nondominated solutions; besides the value represents the distance between two neighboring nondominated solutions of the considered solution. Two particular values remain permanent in the achieve set: the lower and the higher objective function values; their CD possess a significant value and they are also the boundaries. The algebraic form of the CD of the ℎ particle at iteration step is stated as follows: The two particles sent in computation for the th are the − 1 and + 1 particles respectively before and after the considered particle. Therefore, the CD values saved will be compared so that the particles with small CD values are exposed to be replaced for another nondominated solutions. The nondominated solutions altogether which are kept in the external archive set are eligible to be a global best guide (each one). Moreover, a random selection is practicable for all the swarm at each iteration. However, a single steady global best guide occur an efficient convergency. On the other hand, the diversification is jeopardized in all the process. As regards of the approximated Pareto front, all the nondominated solutions recorded in the archive set might be selected to serve as a guide. From the start of the program, a set of particles is run in the PSO variant, then the nondominated solutions are stored in the archive set. Then, at each iteration a nondominated solution is selected as the global best guide randomly. Obviously any nondominated solution can play the role of global best guide. In process of finding the approximated Pareto Front, solutions will converge towards the nondominated solution. According to the parameters' set, the application can guarantee a compelling convergence. Nevertheless the distribution over the archive set remains debatable as in each step all the offspring relay on a random global best guide. In term of perusing the process till the end, many nondominated solutions may gather partially. Therefore, the display of Pareto front could appeal with scattered dense part of nondominated solutions.
Based on the late information, it is suitable to add another tool to hold the diversification. For the current circumstances, a particle with the highest CD value represents the global best guide for the swarm at each iteration. Although the satisfactory was not accurate enough; an additional tool is subsequently added, the DMOPSO-A. This application reinforces the intensification by allowing each particle to possess its own best guide , at each iteration along the entire process. The DMOPSO-A withstands both the intensification and the diversification; by the end the nondominated solutions dwell in a set with flawless diversity characteristic.

Adaptation of genetic operators
During this last decade, operation researchers focus on adapting native evolutionary operators to particle swarm approaches such as mutation [7,27] and elitism [20].

Mutation operator
The mutation operator consists on applying a mutation (a genetic modification) of a particle to improve its position. Applied to a given particle chosen randomly, the mutation operator changes its coordinates with a certain probability. An uniform mutation operator selects with a certain probability a particle to be muted. Then, with another given probability, it chooses a dimension (one of the particle coordinates in the decision space) and apply modification of that position with respect of the search domain space.
In this work, we select the mutation operator [34] which applies the mutation of the personal best guide of the particle instead of the particle itself. The mutation operator is defined by where , is the ℎ coordinate of the personal best guide , of the ℎ particle at iteration step , is a beta distributed random number with parameters , ≤ 1. , is the opposite of , value, i.e., if , = 1, then , = 0 and vise versa. The parameter ′ , is defined by where , is an uniform generated random number named the strategy parameter, means that a normal distributed random number is generated for each dimension coordinate; = 1 √ 2 with is the population size. Note that after the mutation operator applied, comparisons between the personal best guide , and its corresponding muted particle , are made to determine if the new muted particle , dominates , or not.
, is replaced only if its is dominated by its muted particle , .

Crossover operator
In addition to the mutation operator and the elitist policy, the DMOPSO-A approach integrates a crossover operator. This approach, introduced in [34] and named DE-PSO, is formulated as a Differential Evolution algorithm. A target particle is fixed, then a new candidate particle is created by crossing other particles randomly selected in the swarm.
In this work, the crossover algorithm uses three particles selected such that they are pairwise different and all different to the target particle. If denotes the target particle and 1 , 2 , 3 the selected particles for crossover, the crossover operator is defined bỹ︀ where ∼ (0, 1), 1 , , 2, and 3, are the ℎ coordinates of the particles 1 , 2 , 3 at iteration step respectively; and̃︀ , is the ℎ coordinate of the candidate (crossed) particlẽ︀ and 1, is the opposite bit of 1 , value. If the new created particlẽ︀ dominates the target particle , the candidate particle replaces the target one.

Adaptation of the DMOPSO-A for the multiobjective CVRP
Let us describe now the formulation of DMOPSO-A algorithm for the multiobjective CVRP. Scheduling vehicle routing is one of the key issues for a successful application of the DMOPSO-A algorithm to the multiobjective CVRP. The goal is how finding a suitable mapping from problem solution to DMOPSO-A particle. The following modelization has been presented in [4] for the Capacitated Vehicle Routing Problem.
Let us consider a problem of which customers must be visited by vehicles. We can generate a search space with particles positions vectors so that each particle length has × components. Each particle will contain sections, and each section will be formed by discrete points with values 0 or 1. If the value is 1 that means that the corresponding customer is served by the considered vehicle and if the value is 0 this will mean that the vehicle selected does not serve the corresponding customer. An example of this modelization can be set as in Figure 1.
In this example, eight customers must be served by two vehicles. And a possible vehicle routing problem instance (customer, vehicle) is generated. For the decoding process, each particle is represented as a matrix with two lines. For a VRP problem when customers must be served by vehicles, the first line of the matrix is composed by a vector = ( 1 , 2 , . . . , × ) with entries integer values arranged from 1 to × , i.e., = ( = 1, . . . , × ). The route construction is formed as follows: represents a tensor where = 1, . . . , is the candidate particle, the vehicle and the customer. -= 1 if the particle the ℎ vehicle serves the ℎ customer and 0 otherwise. -The vehicle number is computed from the position by using the following equation where ⌊ ⌋ is the integer portion of . -The customer number is extracted using equation -Finally, the scheduled route is obtained by the value contained in the second line of the column . If that value is 1, then the customer is served by the vehicle .
For handling the CVRP constraints, each customer has to be served exactly once by one vehicle, the total vehicle length must not exceed the maximal distance constraint and the total discharge of any route must not exceed the capacity of the vehicle. To check such constraints, the matrix is used. For a given particle 0 , each column (customer service) is examined.
-If there exists only one line such that 0 = 1, then the column (customer) respects the serving constraints.
-If there exists at least two indices (lines) 1 and 2 such that 1 0 = 2 0 = 1 (the customer is served by at least two vehicles), a line 0 is randomly chosen from those with values 1 and its value is maintained while all the others indices will be set to zero. -If the total length in the route exceeds the limited value or the total successive discharge of each direct route exceeds the capacity of the vehicle, the solution is infeasible. For infeasible solutions, carry out the DMOPSO-A operation until the solutions become feasible.

Numerical experiments
In this section, computational results are conducted in order to discuss the performance of the proposed algorithm. Some benchmarks from literature are selected and obtained results are compared to those proposed in literature. All the tests are carried on Matlab R2011b programs under Linux in a PC with AMD E1-1200 APU -2 GB RAM.

Dell'Amico benchmark
The first computational experiment is conducted on the Vehicle Routing problem with Simultaneous Pick-up and Delivery (VRP-SPD) which is a particular case of the Capacitated Vehicle Routing Problem where the vehicle capacity must be revalued after each customer service. Each customer is affected to a demand to be delivered and a given good to be picked. This is a single objective optimization problem where the only need is to minimize the routes cost of the fleet. For instance, the benchmark data set taken into account is recorded in the work of Dell'Amico et al. [11]. This problem comprises five classes of VRP-SPD instances. Each class composed of 20 or 40 customers consists of instances. Class 1 consists of 12 instances, while Class 2S, 2C, 3S, and 3C consist of 18 instances. The following problem parameters are affected in the proposed DMOPSO-A method: fixed cost per vehicle, = 0; variable cost per distance unit, = 1; maximal service duration = ∞; and the number of vehicles available is the number of vehicle in the optimal known solution. The computation test is run during 100 iterations, based on a mutation rate equal to 0.6, a crossover rate equal to 0.39, 1 = 0.75, 2 = 1.2, 3 = 2.6. This problem has been considered in [1] where authors resolve it using a simple PSO algorithm. Their results are using here for comparison.
The comparison of the best solution among 100 iterations of the DMOPSO-A is presented is Table 1. The average total cost of the instances in each class is compared to the one found by Dell'Amico et al. [11]. It has been mentioned that Dell'Amico et al. [11] provided the optimal solution of 75 out of 84 instances and the upper bound represents the optimal solution found. Hence, results in Table 1 demonstrate that DMOPSO-A approach can be able to produce high quality solutions that are very close to the best solution found. We can also remark that only for the 3S class, The Jin approach [1] produces a best solution compared to the proposed approach. In addition, the solutions can be reached in a few computational iterations and a short computation time.

Dethloff test problem
Next, we consider another computational experiment that is based of the benchmark problem presented by Dethloff [12]. This is composed by four data sets named SCA3, SCA8, CON3 and CON8, which consist of 10 instances of a problem of 50 customers. Each one is different of another one by specific characteristics: SCA data are composed by customers distributed uniformly in the service area while the CON problems are generated with half of costumers located uniformly in the service region and the other half are situated in given part of the service affectation. The number after SCA or CON indicated the parameter for determining vehicle capacity.
This VRPSPD is an multiobjective optimization problem where the goal is to minimize the total traveled distance by maximizing the capacity constraint of the vehicle. Therefore, the problem parameters are set as follows: the cost per vehicle is fixed, and = 0; the variable cost per distance unit is = 1; there is no constraint on service duration, so = ∞; and the number of available vehicles is equal to the number of available vehicle in the best known solution.
To compare our approach to those found in the literature, we collect the best solution among 10 iterations based on a mutation rate equal to 0.5, a crossover rate equal to 0.25, 1 = 0.58, 2 = 1.5, 3 = 2.4. The found solutions are compared to those obtained by The Jin [1] using PSO, Tang and Galvao [32] using a Tabu search approach, Bianchessi and Righini [3] by using an heuristic algorithms and the reference solutions in [12]. The overage results of the numerical results over the 10 iterations of each data set are presented in Table 2. The presented results show that the proposed DMOPSO-A performs best results compared the Table 2. Comparison of the 10 iterations overage solution optimal solutions found by DMOPSO-A and others published methods on the Dethloff benchmark.

Set
Average total cost Dethloff [12] Tang [32] Bianchessi [3]   others approaches considered for comparison. Moreover, the proposed MOPSO based approach achieves best results in less computational time. Once more, the considered test problem prove that the proposed approach is competitive and improve solution for this kind of benchmarks.

Solomon's problem
Next, we consider the Solomon's problem which consists of a vehicle routing problem with time windows. The problem is defined as a bi-objective vehicle routing problem with route balancing with 25, 50 and 100 customers. These Solomon's instances are generally used as references in operation research. These are devised into three categories [30,31,39]: category C plus a number of instances where customers are distributed in a cluster, category R plus a number of instances where customers are randomly distributed, and category RC plus a number of instances where the repartition of customers is handled so that distribution is between random and clusters. This problem is solved in [14] using goal programming and genetic algorithm. Their obtained results are used in the current work for comparison. In this example, only the last category with 100 customers is taken     into account. The computation test is run during 10 iterations, based on a population size settled to be equal to 150, a mutation rate equal to 0.5, a crossover rate equal to 0.35, 1 = 0.75, 2 = 1.45, 3 = 1.85.  Experimental results are reported in Table 3 where the first column presents the instance name of the considered Solomon's benchmark [30], the second and third columns point to the best solution (number of vehicles and distance cost) known from literature, the fourth and fifth columns reproduces the average solutions from Ghoseiri's results published in [14]. Finally, the computed best optimal solutions by the current DMOPSO-A approach are listed in the sixth and seventh columns while the eighth and ninth columns represent the average solutions for the 10 runs. The best solution known from literature is taken according to those considered in [14]. Table 8. Two examples of bus orientation found in the nondominated solutions.  Table 8.
These results show that the use of the proposed approach to solve this biobjective optimization problem for the vehicle routing problem based on the Solomon's benchmark is more interesting and competitive compared to results in the literature. This approach lead us to get interesting route balance optimal solutions without considerably affecting the optimal distance cost. We can therefore conclude that the proposed DMOPSO-A approach is adapted to solve VRP with simultaneous pickup and deliveries or with time windows.

Urban school bus routing with route balancing
Finally, we consider an urban school bus routing problem with route balancing. The problem consists of affecting a set of buses to serve a given number of stations to deposit and/or retrieve schoolboys and schoolgirls. In the real-life process, the number of students to be picked or deposed at a given station is stochastic and cannot be fixed a priori. But, in this example, for reasons of simplicity and as preliminary example, this number is considered to be a priori known. This problem is one of the famous applications of the Vehicle Routing Problem with Simultaneous Pick-up and Delivery (VRP-SPD). Because pupils must not arrive late to their school/home, this problem is also one of the applications of the vehicle routing problem with time windows. To take into account the time windows constraint, we use here the route balancing objective. This will lead vehicles to minimize routes circuits and minimize students waiting time in the station and inside the bus. Each bus must return to the depot after the end of service. In this current example, we do not consider the problem of stochastic route traffic which may change the route travel time.
In this example, 20 stations have to be served by at most 10 school buses. To facilitate the modeling, the geographic distribution of the stations must be known. The stations distribution is depicted in Figure 2.
For each station , there are assigned the longitude and latitude localisation as a two-dimensional orthonormal positions, a number of students to be deposed, a number of students to be picked and a service time duration in minutes. The initial delivery number of students are considered to pick-up into the bus from the depot (bus terminal). These are summarized in Tables 4 and 5.
The fleet of buses is considered to be inhomogeneous. So, each vehicle has its own maximal transport capacity and its own maximal service duration. These characteristics are summarized in Table 6.
For comparison, simulations are made for the proposed DMOPSO-A approach. The swarm population contains 20 particles, the mutation rate is equal to 0.75, the crossover rate is setted to 0.6; 1 = 0.48 2 = 0.9 and 3 = 2.0 The code is run up to 500, 1000 and 1500 iterations. Next, Pareto front results corresponding to 1000 and 1500 iterations are plotted in Figure 3.
The compromise solutions are summarized in the Table 7 where the first column of each class corresponds to the total cost of a route while the second column is the route equity (the route balance). Two examples of nondominated particle solutions reported in Table 8 and prove that each route starts from the depot and finishes on the depot; and each station is visited once by exactly one bus. The majority of the buses have almost the same number of stations to serve. This last remark is in concordance of the route balancing desire. These selected two routes balances are depicted in Figure 4.

Conclusion
The Discrete Multiobjective Particle Swarm approach, based on a probability of particle bit update, local attractors, crowding distance, and genetic operators; has been used in this work to solve vehicle routing problems. The performance of this new approach over all the benchmark problems studied turned out to be satisfactory in the sense that the aim was to determine if the new algorithm was able to obtain results close to those underlined in the literature.
Numerical results are shown for three well-known test examples on vehicle routing problems and compared to results from the literature. The obtained results demonstrate that, in spite of its simplicity, handling and flexibility, the DMOPSO-A is a promising approach to multiobjective optimization because its performance was generally better and encouraging. Further work will consider stochastic vehicle routing problems such as Urban School Bus routing without any a priori knowledge of number of students to be picked and/or deposed at a given bus station and by taking into account stochastic traffic flow.