Abstract:
To solve many-objective routing problems, this paper generalizes the Multi-Objective Ant Colony System (MOACS) algorithm, a well-known Multi-Objective Ant Colony Optimization (MOACO) metaheuristic proposed in 2003. This Generalized MOACS algorithm is used to solve a Split-Delivery/Mixed-Fleet Vehicle Routing Problem (SD/MF-VRP) under different constraints, resulting from the mathematical modeling of a logistic problem: the distribution of motorcycles by a Paraguayan factory, considering several objective functions as: (1) total distribution cost, (2) total traveled distance, (3) total traveled time, and (4) unsatisfied demand. Experimental results using the proposed algorithm in weekly operations of the motorcycle factory prove the advantages of using the proposed algorithm, facilitating the work of the logistic planner, reducing the distribution cost and minimizing the time needed to satisfy customers.
Abstract in Spanish:
Para resolver problemas de enrutamiento con muchos objetivo simultáneos, este trabajo generaliza un algoritmo multi-objetivo por Colonia de Hormigas (MOACS), propuesto en 2003. Este algoritmo “MOACS Generalizado” se utiliza para resolver el Problema de Ruteo con Entregas Parciales utilizando una flota heterogénea de vehículos, considerando restricciones no habituales que resultan del modelado matemático de un problema logístico: la distribución de motocicletas por parte de una fábrica en Paraguay, considerando varias funciones objetivo como: (1) el costo total de distribución, (2) la distancia total recorrida, (3) el tiempo total requerido, y (4) la demanda insatisfecha.
]]> Resultados experimentales utilizando el algoritmo propuesto en las operaciones semanales de la fábrica de motocicletas demuestran las ventajas de utilizar el algoritmo propuesto, lo que facilita el trabajo del planificador de logística de la fábrica, reduciendo el costo de distribución y el tiempo necesario para satisfacer a los clientes.Key Words - Multiobjective Optimization; Many-Objective Optimization, Vehicle Routing Problem; Ant Colony Optimization; Distribution.
Keywords in Spanish: Optimización Multiobjetivo; Optimización de Muchos Objetivos, Problema de Ruteo de Vehículos; Optimización por Colonias de Hormigas; Distribución.
Submitted: 2014-11-11. Revised: 2015-04-07. Accepted: 2015-07-12.
1. INTRODUCTION
The Vehicle Routing Problem (VRP) was initially proposed by Dantzig and Ramser [1] in 1959. This problem has been extensively studied due to its enormous importance and its diverse applications in a wide range of practical logistical problems, garbage collection, merchandise distribution, and many other problems that require determining efficient strategies with the objective of minimizing operating costs, or improving the efficiency of the various industrial and commercial processes [2] [3] [4] [5] [6] [7].
There are numerous variants of this paradigmatic problem, arising from the various needs and concrete restrictions of the diverse distribution problems found in real world applications, especially in logistic applications. These variations have been extensively studied in the literature (a complete review may be found in Subramanian thesis [8]). Of all the VRP variants, this article will study a combination of two variants; (1) the Split Delivery VRP (SDVRP), and (2) the Mixed Fleet VRP (MFVRP).
The Split Delivery VRP (SDVRP) generally occurs when the amount of merchandise demanded by a client exceeds the vehicle’s capacity. In this case, the client’s request is satisfied through various vehicles in an attempt to reduce costs, minimize travel time, and/or utilize the fewest number of vehicles [9]. The SDVRP was initially studied by Dror and Trudeau [10], who defined the problem and demonstrated the potential costs savings associated with split deliveries. Archetti, Savelsbergh and Speranza analyze in [9] the maximum achievable cost savings that can be accomplished when utilizing split deliveries. They would later present in [11] a computational study demonstrating that said cost savings depends on the characteristics of each particular case. Some studies demonstrate practical applications of the SDVRP to solve diverse problems ranging from food distribution in a farm [12], helicopter flight schedules [13], or the aforementioned, garbage collection problem [9]. Archetti, Savelsbergh y Speranza in [9] propose a Tabu search algorithm to solve the SDVRP, while in Sui, Tang, and Liu’s article [14], it is demonstrated that the split delivery problem can be efficiently solved through an Ant Colony Optimization (ACO) algorithm, which will also be used in this work .
The Mixed Fleet VRP (MFVRP), implies vehicles with different capabilities (or heterogeneous capabilities), with known fixed and variable costs related to each vehicle in a fleet that must serve a series of consumers with known demands. In [15], Golden, Assad, Levy and Gheysens describe a series of effective heuristic procedures for the problem of routing with a heterogeneous fleet, with the objective of determining the optimal truck fleet size and its capabilities, minimizing a cost function. The authors, Subramanian, Penna, Uchoa, and Ochi [16] studied the optimal composition of a fleet of vehicles through a hybrid algorithm, as well as determining the routes that would minimize travel expenses. Similarly, Salhi and Rand, in [17], and Taillard in [18], also attempt to find the ideal composition for a fleet of vehicles by solving the MFVRP. The authors of [19], Wassan and Osman, developed new Tabu Search (TS) variants in order to solve the heterogeneous fleet problem. At the same time, the article by Chen and Ching [20] suggests the alternative of employing an Ant Colony Optimization algorithm in order to solve the heterogeneous fleet routing problem, proving that ACO is a competitive algorithm for this VRP variant; another factor influencing its adoption for this work.
The present study is set upon solving a merchandise distribution problem, in which m establishments must be supplied from a central warehouse, utilizing a heterogeneous fleet consisting of k vehicles always leaving from the same central warehouse, where orders can be fulfilled in more than just one trip (split delivery). This particular problem is called SD/MF-VRP, and it has already been studied in the specialized literature, for instance, by Belfiore and Yoshizaki [21] and [22], in a mono-objective context; i.e., optimizing a single objective function at a time. However, this work proposes to solve for the first time, an SD/MF-VRP variant, taking into consideration the simultaneous optimization of several objective functions in a purely multi-objective context, where no objective function is necessarily more important than the others, facing a variant of the SD/MF-VRP problem with a concrete utility in a motorcycle factory in Paraguay, considering specific third-world restriction, for example, that not all vehicles may transit on any road due to the condition of the various routes, or because of limitations related to weight and/or height.
In order to solve the aforementioned SD/MF-VRP multi-objective problem, we propose utilizing one of the most successful metaheuristics to solve routing problems [23], the Ant Colony Optimization (ACO), which consists of a technique that simulates the indirect communication utilized by ants to establish the shortest path from their nest to its food source and back. Said metaheuristic was first introduced with the name of Ant System (AS), by Dorigo, Maniezzo, and Colomi [24], contributing to a great extent in the development of new computational techniques in the optimization area. Subsequently, in [25], Dorigo and Gambardella presented the Ant Colony System (ACS), applying it to the Traveling Salesman Problem (TSP), proving that this technique turns out to be very competitive and even exceeds other bio-inspired algorithms, such as generic algorithms and other evolutionary variants [23].
]]> When considering a multi-objective VRP approach, it is worth mentioning the proposal put forth by Baran and Schaerer [26], inspired in the Multiple Ant Colony System (MACS) for the Vehicular Routing Problem with Time Windows (VRPTW). Indeed, while the MACS algorithm utilized two ant colonies to minimize two different objective functions (number of vehicles and total travel distance), Baran and Schaerer presented in [26] a Multi-Objective Ant Colony System (MOACS), which utilized a single colony and considers a differentiated “visibility” scheme for each ant in order to improve the exploration by way of making each ant focused its search in different areas of the definition of the problem domain. This first proposal by Baran and Schaerer [26] considered only 2 objectives, solving the Biojective-TSP; therefore, in what follows the MOACS will be generalized to a larger number of objectives.Other articles propose the utilization of ACO to solve variations of the VRP. Lezcano, Pinto, and Baran [27], present experimental results on a test bed of three different families of problems: the (1) Traveling Salesman Problem (TSP), (2) Quadratic Assignment Problem (QAP), and (3) Vehicle Routing Problem with Time Windows (VRPTW), proposing a new approach known as Distributed Team Ant Algorithms for the solution of optimization problems considering various objectives. On the other hand, Hermosilla and Baran [23] performed a comparison between an Ant Colony System, and an evolutionary strategy for the resolution of multi-objective vehicle routing problem with time windows, evidencing that the Ant Colony System has in general a better performance when compared to other evolutionary strategies, especially for large problems.
It is worth mentioning that there are many other proposals of Multi-Objective Ant Colony Optimization (MOACO), such as [27], [28], [29], [30], [31] and [32], but this particular work will be based on the MOACS [26], given its clearly competitive characteristics when compared to other known alternatives [33].
Finally, it is worth mentioning that this work proposes for the first time to solve the SD/MF-VRP considering four simultaneous objective functions: the (1) number of vehicles, (2) total travel time, (3) total delivery time, and (4) unsatisfied demand, making the problem especially complex, given that it is already known that increasing the number of objective functions will also considerably increase the problem complexity [34]. In particular, it is worth mentioning that mainly because of this increasing complexity, when the number of objectives is larger than three (as it is the case with the problem at hand), the multi-objective optimization problem even receives a special name: Many-objective Optimization Problem [34]. This special consideration entails the need of rethinking the MOACS algorithm to be implemented when 4 or more objectives are considered, giving rise to a need of generalizing the proposed algorithm [26], as presented in Section IV.
In conclusion, this work studies the resolution of a variant to the SD/MF-VRP problem, emerging from a practical, real case in a motorcycle factory, taking into account four simultaneous objective functions (i.e., considering a many-objective problem), proposing a generalized version of the MOACS algorithm which is able to deal with any number of objective functions (even 4 or more objectives). The rest of the article is organized in the following manner: next section very briefly presents the multi-objective optimization problem, whereas the Mathematical Model of the problem to be explored is described in section 3. Section 4 presents the optimization technique using an ant colony optimization, leaving the specific algorithm proposed for this work to Section 5. The application of the algorithm to the distribution of motorcycles in a factory in Paraguay is found in section 6, while the experimental results are found in section 7, leaving conclusions to section 8.
2. MULTIOBJECTIVE OPTIMIZATION
As its name suggests, in a multiobjective optimization problem, each solution is compared considering more than one objective function. Each of these objective functions must be minimized or maximized. One general multiobjective optimization problem includes a series of n decision variable (n ≥ 1), u objective functions (u ≥ 2), and l restrictions (l ≥ 0). The objective functions and restrictions are functions of the decision variables. This can be expressed as follows:
where 
is the decision vector (independent variable), while 
is the objective vector.
Depending on the type of problem, optimizing can mean minimizing or maximizing. Moreover, some objective functions may be minimized while the rest of the objective functions are maximized. Without loss of generality, a pure minimization context will be assumed in the rest of this article.
In general, there is no one unique optimal solution for a multiobjective problem considering various objectives; rather, there may be a set of compromise solutions that cannot be considered neither better nor worse than the other solutions in this set when all objectives are considered simultaneously. Because of this, in a multiobjective optimization context, the concept of dominance is used. This way, it is said that one solution dominates another if it is not worse in any objective function and it is strictly better in at least one objective. The set of all feasible solutions (meeting all restrictions) not dominated by other feasible solutions, make up the optimal solution set of a multiobjective problem and it is known as a Pareto set. The projection of a Pareto set to its corresponding objective space is known as the Pareto front [26].
3. MATHEMATICAL MODEL
]]> Let,
be the set of branches, where s0 stands for the warehouse and si represents branch i
the set of total capacities of each of the v vehicles available, where Qk is the capacity of vehicle k
Z = {zkij} indicates if vehicle k can transit from si to sj (in which case zkij = 1) or it cannot (zkij = 0)
dij the distance from si to sj. Only the symmetric case is considered in what follows, i.e., dij = dji
tij travel time from branch si to sj. This work assumes tij = tji
qi demand of branch si for 1≤ i ≤ m. It is assumed q0 = 0
C (k, i, j) the cost of using the vehicle k to go from si to sj
a solution to the problem, whose elements are the matrices 
; i.e. 
explained later.
Minimize the vector F(ψ) = [F1(ψ), F2(ψ), F3(ψ), F4(ψ)]T
where F1(ψ), F2(ψ), F3(ψ) and F4(ψ) represent the four objective functions that are next defined.
Total Cost F1(ψ)
(1)
Total travel time F2(ψ)
(2)
Total distance traveled F3(ψ)
]]> 
(3)
Unsatisfied demand F4(ψ)
(4)
The problem is subject to the following restrictions:
Restriction 1
The capacity of a vehicle cannot be exceeded:
(5)
indicates the number of motorcycles carried by vehicle k to branch i. Note that 
, i.e., the number of motorcycles carried by vehicle k to branch si cannot be greater than the order 
made by branch i.
Restriction 2
Vehicles can only transit on permitted roads.
xkij = 1 only if zkij = 1 (6)
where
A solution 
is represented by two matrices with a dimension of 
, where Nmáx indicates the maximum number of branches that may be visited in a single trip.
Represents the i-th branch visited by vehicle k.
Represents the total number of motorcycles that a vehicle k carries to the branch 
.
4. ANT COLONY OPTIMIZATION
]]> Ant Colony Optimization (ACO) is one of the most recognized meta-heuristics for the studied problem given its already known efficiency in solving routing problems [23]. The optimization procedure used by the Ant Colony Optimization is inspired in the behavior of biological ant colonies in order to solve optimization problems utilizing traces of pheromones. This behavior used by biological ants has inspired this metaheuristic to solve optimization problems [30] utilizing colonies of artificial ants, meaning, simple computational agents working in a cooperative manner and indirectly communicating through artificial traces of pheromones.The ACO meta-heuristics are essentially constructive algorithms: at each iteration of the algorithm, each ant (or agent) builds a solution to the problem touring a graph of construction [26]. Every edge of the graph, representing every step an ant could take, has two types of information associated with it to guide an ant’s movement:
Heuristic Information, that measures the preference of moving from node i to another node j of the graph, meaning, touring the link (i, j). This heuristic information is denoted as ηij and it is known as visibility.
Information about traces of artificial pheromones, that measures the “learned desirability” of a movement from a node i to another node j, deposited by the ants. This information is modified during the execution of the algorithm depending on the solutions found by the ants; therefore, they store the learned information during the resolution process. It is denoted as τij and it is known as pheromone.
Different versions of Multi Objective Ant Colony Optimization (MOACO), differ from each other mainly in three points: the (1) pheromone traces, (2) solutions to be rewarded, and (3) determination of heuristic factors.
Pheromone traces. The quantity of pheromone in a component represents the colony’s experience related to choosing that component. When there is only one objective function, this experience is defined using this sole objective. However, when there are several objectives, at least two strategies can be considered. One strategy consists in a single pheromone structure, in which case, the quantity of pheromone left by ants is defined according to an aggregation of the different objectives considered when solving the problem at hand (each aggregation option, as a weighted sum, can give rise to a different algorithm). The second strategy tries to consider various pheromones structures, one for each objective function.
Solutions to be rewarded. When the pheromone traces are updated based on what was learned by the algorithm, it must be decided which of the constructed solutions will leave pheromones, increasing the probability of using the arcs of the solution reached once again. One possibility is to reward solutions that find the best values for each criterion within the current generation. Another possibility is to reward each no-dominated solution within the current generation. In this case, every solution in a Pareto set can be rewarded, or perhaps only the new solutions entering the Pareto set in the current generation.
Determination of the heuristic factors. When solutions are constructed, one candidate is chosen according to a probability depending on a pheromone factor and a heuristic factor (known as visibility) which measures the preference of moving from one node to another. To define the heuristic factor, at least two strategies could be considered. The first assumes an aggregation of the different objectives within a single heuristic information. A second strategy requires considering each objective separately. In this case, there is usually a different colony for each objective. However, competitive results were also reported when a single colony uses a different visibility for each objective function [26].
5. MOACS GENERALIZATION
]]> This work draws inspiration from the approach proposed by Baran and Schaerer in 2003 in order to solve the bi-objective TSP [26], which consists in a MOACO algorithm with a single ant colony, a single pheromone structure and different adaptive visibility for each objective function, known as MOACS (Multi-objective Ant Colony System), considering the promising results from experimental tests that this algorithm achieves when compared with other alternatives when solving a multiobjective problem [33].The MOACS utilizes a single ant colony in order to simultaneously minimize every objective function, considering all objectives as equally important; therefore, a whole Pareto-optimal solution set can be found in only one run of this meta-heuristic. Every objective share the same pheromone trace, represented in a single matrix τ = { τij }. The MOACS algorithm’s central idea applied to the motorcycle factory problem explored in this work is to build only feasible solutions utilizing as many vehicles as needed, given that the minimization of the number of vehicles is not a desired objective.
At each generation (or iteration) of the proposed algorithm, every ant h (of a set of H ants) constructs a feasible solution, beginning its tour in a warehouse (s0) and successively choosing a following node (or warehouse to return) θj, out of the set of feasible nodes not yet visited 
where the subindex i indicates that the ant h is in the node θi. At each node θi, the set is calculated in order to discover every node (or client) to which merchandise still needs to be taken and are not in violation of any restriction (such as a vehicle size, or that is too heavy to travel through a given road). The set does not include the warehouse until a tour has been finished. Naturally, when a vehicle cannot visit any more nodes, for instance because it is fully loaded, it must return to the warehouse and another vehicle starts choosing nodes in the same manner. This process is repeated until every client has been served and therefore, a solution has been found.
An ant h choses to move from one node θi to another node θj using both: the heuristic information and the pheromone traces. The heuristic information is given by the visibility 
, while the pheromone information is given by τij. To choose the next node θj to be visited by an ant h which is found in θi, each node or branch is assigned a probability pij giving in (7), and, in this manner, the next branch θj is chosen from 
at random using that probability.
It is relevant to remember that the excellent results obtained by the MOACS are due in great part to the use of different visibilities, one for each objective function, and a mechanism for each ant h to prioritize its search in different parts of the domain of feasible solutions, which is why this characteristic originally given for only 2 objectives must be generalized here for a number u of objective functions. In consequence, the probability pij can be generalized to u objectives in the following manner:
(7)
where,
u is the number of objective functions of the problem, while 
is the visibility related to the objective function u for the link (i, j);
the variables 
define the relative influence among visibilities; and
parameter α defines the relative influence of the pheromone traces.
]]> When generalizing Baran and Schaerer’s original proposal [26] for u objectives (u ≥ 2), each ant utilizes one of the 
different combinations for the visibility, as shown in Table I, i.e. 
, in order to improve the exploration of the search space. Consequently, it is recommended to use 
ants in order to cover all possible variants found in Table I. One example with 
is presented in Section VI, and it is used to solve the studied SD/MF-VRP problem for a motorcycle factory.
The value taken by each variable of relative influence 
varies with the number h of the ant building a solution, achieving in this way the desired characteristic of having each ant h having preferences on different areas of the exploration domain. This way, each ant h potentially aims at different regions of the Pareto front that is being searched.
In consequence, the proposed solution considers generations in which H ants construct solutions, each with a different probability distribution, expecting in this way to improve the exploration of the MOACO.
TABLE I. Possible values for the visibility parameter λ utilized by the different ants in calculating probabilities.
When each ant h of a generation w finds a complete solution 
, this is compared with the already stored solutions in the best known Pareto set approximation to verify if it is a no-dominated solution. If, in effect, it is a new optimal Pareto solution, it is included in a set P and dominated solutions are erased from the current set P of optimal solutions.
After each generation of the proposed algorithm is finalized, an update of the pheromone matrix is carried on with the no-dominated solutions contained in the Pareto_Set, according to the equations
τ’ij = τij + Δτij (8)
]]> τ’’ij = (1 - ρ).τ’ij (9) For this calculation, the objective functions are normalized between a very small value 
and 1; this is, to the interval 
for instance, utilizing the following equation:
(10)
Therefore, to calculate the value of Δτij the following equation may be used:
(11)
where 
is a weight corresponding to the relative importance that can be attributed to each objective function 
and that can be defined at the discretion of the decision maker or at random each time it is used. In the actual case presented in section VI, the selection was done by the chief of logistics of the motorcycle factory that implemented the presented generalized MOACS.
The objective of the proposed algorithm is to only build feasible solutions, utilizing the strictly necessary vehicles, until every client has been served or until the capacity of the vehicle fleet has reached a saturation point. A finite number of vehicles is considered for distribution, thus, it will not always be possible serving every branch. Indeed, the algorithm can also end when the capacity of all available vehicles has been met, even if there are still branches to visit. The unsatisfied demand from the branches left out will be considered as another objective function to be minimized, and, in the particular case of the motorcycle enterprise, could be satisfied the following day.
Figure 1 shows the general procedure for the MOACS algorithm implemented in this work, remembering that at each generation w, H ants construct solutions 
, calculating their respective objective functions:
.
is not worse than another solution 
(line 3) in any objective function and also, the solution is strictly better than solution in at least one objective (line 4), then it is said that the solution is dominated by . The dominated is marked and discarded (line 5), while the no-dominated solutions are stored in a set P of no-dominated solutions. This procedure is used in the algorithm in two instances: to find the no-dominated solutions at each generation, and then, to update the Pareto_Set with these no-dominated solutions. 
Figure 1. Optimization Algorithm implemented
6. APPLICATION TO MOTORCYCLE DISTRIBUTION
The method proposed in section V has been applied to an actual case utilizing real and historical data concerning a motorcycle manufacturing enterprise in Paraguay [35].
The motorcycle factory studied did not apply a formal method in the planning of its vehicle routes when distributing their products to the different branches, thus, this job was naturally tedious for the employee in charge of logistics, who was satisfied enough with being able to automatize the procedure as much as possible. The logistic area of the factory worked in the following way: the department in charge of logistics within the business collected weekly orders from their internal clients (branches) and continuously made empirical decisions without a mathematical model that would allow them to neither quantify their true costs nor take decisions that would allow the enterprise to optimize their distribution. In consequence, this work mathematically models the logistical problem with the distribution of motorcycles and proposes the utilization of the Generalized MOACS algorithm presented in the previous section to solve the presented mathematical model.
Figure 2. Procedure to find no-dominated solutions
This particular company has at hand a central warehouse established at the factory site and m = 58 branches that must be supplied to, from the central warehouse. The company, besides, has a fleet of v = 5 heterogeneous vehicles destined exclusively to the distribution of supplies, always leaving from (and returning to) the central warehouse (the factory).
After talking with logistics experts, it was decided to simultaneously minimize the following four objective functions:
: total merchandise distribution cost,
: total traveled distance, 
: total travel time, and
: unsatisfied demand in a given day.
To calculate the probability that a branch sj will be visited from a branch si, equation (7) found in Section V is used. In said equation, 
represents the set of feasible nodes or branches that have yet to be visited (excluding the warehouse) and not in violation of any restriction to the problem. In this probability formula, as many visibilities as there are objective functions considered in the formula are used; however, for the actual application to the motorcycle factory it was decided to use only three visibilities related to the objective functions of cost, distance, and time. This is because, for the particular case of the business being studied, throughout the length of a work week, demand can always be satisfied by hiring third party vehicles if it is required; therefore, the visibility corresponding to the unsatisfied demand was not considered.
The three visibilities used to calculate the probability are defined as:
refers to the visibility related to the cost cij
refers to the visibility related to time tij
refers to the visibility related to distance dij
Finally, the formula used to calculate the probability of visiting a branch sj from a branch si, is the following:
where,
parameter α defines the relative influence of pheromone. For this work, α=1 is chosen,
variables 
define the relative influence among visibilities.
For the tests performed, the variables assume values of 0, 1 y 2, that is, 
, meaning that . Table II clearly shows all 27 combinations of the values for the parameter used in the application that was implemented.
The value assumed by each variable of relative influence varies with the number h of the ant building a solution, as it can be observed in Table II. It is worth noting that when = 0, the corresponding visibility has no effect in the probability value.
In consequence, the solution proposed to the motorcycle factory considers generations in which 27 ants construct solutions, each one with a different probability distribution, hoping that way to improve the MOACO exploration.
]]>TABLE II. Parameter λ values used in the final implementation for the motorcycle factory in Paraguay
For the updating of the pheromone matrix τ, a Δτij has been considered, uniquely associated with the distribution cost, at the explicit request of the logistics expert, due to the fact that the motorcycle company ensures that cost minimization is their priority, denoted as . Therefore:
(13)
It is worth noting that the choice of the method to calculate Δτij can affect the quality of the calculated results, prioritizing some objectives, but under no circumstances it deletes from the algorithm the ability to find a set of different multiobjective solutions that make up the Pareto set.
The entry parameters considered in the application of the algorithm are detailed next:
G (V, A) represents a graph; that is, the set of nodes (or branches), linked by edges (or roads);
α defines the relative influence of the pheromone traces;
]]> parameters define the relative influence between visibilities; where 
represents the relative influence of the cost visibility. Analogically,
and 
represent the relative influence of the total travel time and total distance travelled visibilities respectively; ρ represents the evaporation rate ;
D = {dij} is a matrix with a dimension of (m+1) x (m+1) that includes the central warehouse and the m company branches. Each element dij from D, indicates the distance between a branch si and another branch sj, in km. It is considered that dij = dji. This matrix is used to calculate the objective function , representing the total distance traveled in a solution ; also, it is used for calculating the relative visibility to the distance traveled 
.
T = {tij} is a matrix with a dimension of (m+1) x (m+1) that includes the central warehouse and the m company branches. Each element tij represents the travel time between a branch si and another branch sj, in hours. It is considered that tij = tji. The matrix T is used for calculating the objective function that represents total travel time as well as for calculating the relative visibility to time 
.
R = {rij} is a restriction matrix with a dimension of (m+1) x (m+1) given for each vehicle k that has circulation restrictions (a single vehicle for this particular factory), that includes the central warehouse and the m company branches considered for this work. Each element rij on the R matrix indicates whether a vehicle k with restrictions, can tour (or not) the link (i,j), where,
.
- ]]>
qi represents each branch’s si demand;
Cost_veh is a column vector, containing the parameter cost/km for each vehicle k, where 
. The vehicles cost/km parameter is multiplied by the distance traveled by each vehicle (in km), to calculate the objective function , representing a solution’s total cost.
To calculate the visibility relative to the cost, 
, matrices Ck = {ckij} must be considered, where . Every matrix Ck must include the traveling costs between a branch si and another branch sj, when a vehicle k is used. Every element ckij would be attained by multiplying the parameter cost/km of the corresponding vehicle, by each element dij in the distance matrix D.
7. EXPERIMENTAL RESULTS
Experimental testing was first conducted considering the historical orders and delivery records of the Paraguayan motorcycle company. In such records it could be found data including: order and delivery dates, demand from each branch, quantity of motorcycles delivered to each branch, and which vehicle visited each branch. For the tests next reported, requests from all 58 company branches were considered.
Figure 3 indicates the weekly procedure followed by the motorcycle company when using the software developed to determine the distribution of received orders. The logistic of the company is organized in such a way that orders mostly arrive on Thursdays and Fridays, but unexpected or urgent requests may continue arriving during the first days of the following week. The software developed for the motorcycle factory is executed for each one of these days, having as input the demands of all branches. As a result, the software proposes the planning for the distribution of these orders, to be carried on the following day.
Logically, it could happen that a day’s demand exceeds a fleet of available vehicle’s capacity, considering both, their own vehicles and those of a third party. In this case, unsatisfied demand is added to the new orders received on that day, and this new total demand is reentered the program to plan next day distribution. This process continues until every order received in the week has been delivered.
Tables III, IV and V, present a comparison of the actual distribution performed by the company with a solution of the Pareto set calculated with the proposed method in order to emphasize the advantage associated with using the proposed Generalized MOACS metaheuristic. Comparisons are made in terms of the total cost in Gs. (official currency of Paraguay, where 1 US Dollar ≈ 5000 Gs.), the total traveling time, and the total traveled distance, as well as a comparison of the quantity of days needed for the entire demand to be distributed (last column, called Distribution days). The comparisons do not include unsatisfied demand due to the fact that in every test performed it turned out that it had no effect given that the vehicle fleet is large enough to satisfy the internal customers of the factory, as was already put forth when explaining the simplifications made in the algorithm having this fact in mind. All data necessary to reproduce the results seen in tables III, IV and V, is available in [35].
For the first studied period (see Table III), it was found that the distribution actually made at the factory is part of the Pareto set, which proves the high level of specialization of the logistics expert. However, with the method proposed in this work, at least one lower cost solution is obtained, even though, logically, the objective functions of time and distance are better when utilizing the company’s actual distribution given it is an optimal Pareto solution (see Table III). For the second period considered (see Table IV), one solution found by the Generalized MOACS clearly dominates the distribution programmed by the company specialist. In fact, the calculated solution strongly dominates the human expert solution, i.e. the calculated solution is better in every aspect of Table IV; clearly demonstrating the usefulness of the proposed method.
]]> Figure 3. Weekly procedure for determining the distribution of orders at a Paraguayan factory of motorcycles.
Table V shows how, in the third test carried on, the proposed method achieves a solution which strongly dominates the distribution programmed by the company, given that, once again, it is better in every compared aspect. In fact, for the this final studied period, the proposed algorithm finds a solution that significantly minimizes the three objectives with respect to the company’s actual distribution system. On top of that, the proposed method also enables branches demand to be satisfied in two days less than what was achieved by the company’s logistics department.
TABLE III. FIRST COMPARISON TEST
TABLE IV. SECOND COMPARISON TEST
TABLE V. THIRD COMPARISON TEST
The three experimental tests conducted demonstrate that the proposed Generalized MOACS method achieves good feasible solutions for the company’s distribution, generally better than the actual solutions found with the empirical method presently used by the company’s logistics expert who dedicates a considerable fraction of his working time to organize the motorcycle distribution. Also, it is worth noting that in general, most compared objectives are simultaneously minimizing, i.e. the total cost, total traveled distance, total traveled time and unsatisfied demand may be improved at the same time, keeping in mind that in two out of three conducted tests, algorithm solutions strongly dominate (in all objectives considered) the company’s logistics expert planned distribution system. At the same time, the program found solutions capable of reducing costs in all three conducted tests, a fact that seems relevant for the company. This makes the proposed method highly attractive for the business, considering the owners priority of minimizing costs. Another important advantage of the proposed method is that it considerably reduces the time needed to plan the distribution, when going from a manual programming to a computerized one.
8. CONCLUSIONS AND FUTURE WORK
]]> When trying to solve a Split Delivery/ Mixed Fleet – Vehicle Routing Problem or SD/MF-VRP for a motorcycle distribution company in Paraguay, we found ourselves with a need to consider various unconventional restrictions (such as vehicles that cannot transit on certain roads) and to simultaneously minimize 4 objectives: (1) the total cost, (2) total travel time, (3) total delivery time, and (4) unsatisfied demand, which entails solving a pretty complex practical problem, today known as many-objective optimization problem [34]. In consequence, after carefully analyzing the state of the art, a well-known MOACO was chosen to solve the problem at hand, the MOACS [26] proposed in 2003 to solve a bi-objective TSP problem. In consequence, said algorithm was modified to be able to treat a generic number u of objective functions (u ≥ 2). This generalized version of the MOACS explained in detailed in section V, was then successfully tested at a Paraguayan motorcycle factory, as explained in section VI, proving the viability of the proposal.Experimental results using historical data and later validation at the motorcycle factory that requested the developed solution allow us to confirm that the proposed Generalized MOACS achieves good feasible solutions, similar or even better than the factory’s empirically solutions obtained by its logistic expert. It was also proven in the three conducted tests that if the company had planned their distribution with the proposed algorithm, they could have achieved a significant costs savings in the distribution final cost. In the first and second tests, a cost savings of Gs. 660,000 could have been generated, while for the third test it was calculated that, had the company used the proposed algorithm, the saving would have been in the order of Gs. 3,400,000.
After day to day trials at the motorcycle factory, it was also clear that the proposed method manages to optimize the available vehicle fleet’s utilization, minimizing third party vehicles, increasing in this way the use of their own vehicles.
It is worth noting that the company’s Logistics department does an excellent job in planning distribution routes; nevertheless, being an empirical planning, there is a lot more room for human mistakes to be made, which the proposed algorithm is capable of systematically avoid, reducing time and letting the expert choose from a set of optimal Pareto solutions, what helps the logistic worker to better understand the trade-off among different alternatives. Even the logistic expert said he learned from the possibility of knowing the possible trade-off solutions without the huge work he needed before automation using the developed system.
Considering the results from the conducted tests, it is finally concluded that the proposed method, which programs the distribution of motorcycles according to the MOACO algorithm represented in Figure 1 and complemented in Figure 2, constitutes an excellent alternative for the planning of distribution routes for the company in question, what clearly can be imitated by other companies with similar logistic distribution problems.
In order to continue with the work already under way, the following topics are suggested as future jobs:
Considered that some branches could have their own, small warehouse, which could potentially satisfy orders made in their vicinities. Improve the problem’s approach by considering that motorcycles can be delivered from another branch and not necessarily from the central warehouse.
Add new objective functions and/or restrictions, for example, establishing that all drivers should travel a similar amount of daily or weekly hours.
- ]]>
Consider the motorcycle stock problem. This way, before programming the distribution, there must be an assurance that the ordered models are in fact available in the warehouse.
Consider a factory where motorcycle may have different sizes.
Establish the programming period to be configurable: daily, weekly, or monthly.
Apply the Generalized MOACS algorithm to other many-objective problems with several objectives. Consider the generalization of other evolutionary algorithm, even any known MOACO, using similar techniques to the one presented in this work.
REFERENCES
[1] G. B. Dantzig and J. H. Ramser, “The Truck Dispatching Problem,” Management Science, vol. 6, nº 1, pp. 80-91, October 1959.
[2] P. Toth and D. Vigo, The Vehicle Routing Problem (Vol. 9), Philadelphia: SIAM, 1987, pp. 1-26.
[3] M. Vásquez, “Desarrollo de un Framework para el Problema de Ruteo de Vehículos,” Thesis at Universidad de Chile, Santiago de Chile, 2007 (www.tesis.uchile.cl/handle/2250/102934).
[4] B. Golden, S. Raghavan and E. Wasil, “The Vehicle Routing Problem Latest Advances and New Challenges,” de Operations Research Computer Science Interfaces, New York, Springer, 2008.
[5] G. Laporte, “The vehicle routing problem: An overview of exact and approximate algorithms,” European Journal of Operational Research, vol. 59, Nº 3, pp. 345-358, 1992.
[6] G. Laporte, M. Gendreau, J.-Y. Potvin and F. Semet, “Classical and modern heuristics for the vehicle routing problem,” International Transactions in Operational Research, vol. 7, nº 4-5, pp. 285-300, 2000.
[7] N. Christofides, “The vehicle routing problem,” RAIRO - Operations Research, vol. 10, nº 1, pp. 55-70, 1976.
[8] A. Subramanian, Heuristic, Exact and Hybrid Approaches for Vehicle Routing Problems, Thesis at Fluminense University - Brazil, 2012.
[9] C. Archetti, M. W. P. Savelsbergh and M. G. Speranza, “Worst-Case Analysis for Split Delivery Vehicle Routing Problems,” Transportation Science, vol. 40, pp. 226-234, 2006.
[10] M. Dror and P. Trudeau, “Savings by Split Delivery Routing,” Transportation Science, vol. 23, nº 2, pp. 141-145, Mayo 1989.
]]> [11] C. Archetti, M. Savelsbergh y M. G. Speranza, “To Split or Not to Split: That is the Question,” Transportation Research Part E: Logistics and Transportation Review, vol. 44, pp. 114-123, Enero 2008.[12] P. Mullaseril, M. Dror and J. Leung, “Split-delivery routing in livestock feed distribution,” Journal of the Operational Research Society, vol. 48, pp. 107-116, 1997.
[13] G. Sierksma and G. A. Tijssen, “Routing helicopters for crew exchanges on offshore locations,” Annals of Operations Research, vol. 76, nº 0, pp. 261-286, 1998.
[14] L. Sui, J. Tang y S. Liu, Ant colony optimization algorithm to solve split delivery vehicle routing problem, Yantai, 2008, pp. 997-1001.
[15] B. Golden, A. Assad, L. Levy y F. Gheysens, “The Fleet Size and Mix Vehicle Routing Problem,” Computers & Operations Research, vol. 11, Nº 1, pp. 49-66, 1984.
[16] A. Subramanian, P. Penna, E. Uchoa and L. Ochi, A Hybrid Algorithm for the Fleet Size and Mix Vehicle Routing Problem, Metz, 2011.
[17] S. Salhi and G. Rand, “Incorporating vehicle routing into the vehicle fleet composition problem,” European Journal of Operational Research, vol. 66, Nº 3, pp. 313 - 330, 1993.
[18] É.D. Taillard, “A heuristic column generation method for the heterogeneous fleet VRP,” Lugano, Switzerland, 1996 (http://mistic.heig-vd.ch/taillard/articles.dir/Taillard1999.pdf).
]]> [19] N. Wassan and I. Osman, “Tabu Search Variants for the Mix Fleet Vehicle Routing Problem,” The Journal of the Operational Research Society, vol. 53, Nº 7, pp. 768-782, Julio 2002.[20] C.-H. Chen and C.-J. Ching, “An Ant Colony Optimization Algorithm for the Heterogeneous Fleet Vehicle Routing Problem,” Journal of the Eastern Asia Society for Transportation Studies, vol. 9, pp. 631-643, 2011.
[21] P. Belfiore and H. T. Yoshizaki, “Heuristic methods for the fleet size and mix vehicle routing problem with time windows and split deliveries,” Computers & Industrial Engineering, vol. 64, Nº 2, pp. 583-601, 2013.
[22] P. Belfiore and H. T. Yoshizaki, “Scatter search for a real-life heterogeneous fleet vehicle routing problem with time windows and split deliveries in Brazil,” European Journal of Operational Research, vol. 199, Nº 3, pp. 750-758, 2009.
[23] A. Hermosilla and B. Barán, Comparación de un Sistema de Colonias de Hormigas y una Estrategia Evolutiva para un Problema Multiobjetivo de Ruteo de Vehículos con Ventanas de Tiempo, Arequipa, 2005, pp. 379-388.
[24] M. Dorigo, V. Maniezzo and A. Colorni, “The ant system: An autocatalytic optimizing process. Technical Report 91-016,” Milan, 1991.
[25] M. Dorigo and L. M. Gambardella, “Ant Colony System: A cooperative Learning Approach to the Traveling Salesman Problem,” Evolutionary Computation, IEEE Transactions, vol. 1, pp. 53-66, Abril 1997.
[26] B. Barán and M. Schaerer, “A Multiobjective Ant Colony System for Vehicle Routing Problem with Time Windows,” de Applied Informatics, pp. 97-102, 2003.
[27] C. Lezcano, D. Pinto, and B. Barán, Team Algorithms Based on Ant Colony Optimization–A New Multi-Objective Optimization Approach, En Parallel Problem Solving from Nature–PPSN X. Springer Berlin Heidelberg, 2008. pp. 773-783.
]]> [28] D. Pinto, “Enrutamiento Multicast Multiobjetivo Basado en Colonia de Hormigas,” Thesis at the National University of Asunción, Paraguay, 2005.[29] A. Arteta, B. Barán and D. Pinto, “Routing and wavelength assignment over WDM optical networks. A comparison between MOACOs and classical approaches,” de IFIP/ACM Latin-American Networking Conference - LANC'2007, San José de Costa Rica, 2007.
[30] C. Lezcano, D. Pinto and B. Barán, “Team Algorithms based on Ant Colony Optimization. A new Multi-Objective Optimization approach,” de 10th International Conference on Parallel Problem Solving From Nature - PPSN'2008, Alemania, 2008.
[31] D. Pinto and B. Barán, “Multiobjective Multicast Routing with Ant Colony Optimization Approach,” en Proceedings of the 3rd international IFIP/ACM Latin American conference on Networking. ACM, 2005. p. 11-19.
[32] C. Solnon and K. Ghédira, “Ant colony optimization for multi-objective optimization problems,” 2012 IEEE 24th International Conference, vol. 1, pp. 450-457, 2012.
[33] C. García-Martínez, O. Cordón and F. Herrera, “A taxonomy and an empirical analysis of multiple objective ant colony optimization algorithms for the bi-criteria TSP,” European Journal of Operational Research, vol. 180, Nº 1, pp. 116-148, 2007.
[34] C. von Lücken, B. Barán and C. Brizuela, “Survey on multi-objective evolutionary algorithms for many-objective problems,” Computational Optimization and Applications, pp. 1-50, 2014.
[35] G. Insaurralde and M. Laufer, Tesis de grado: Método formal de planificación de rutas de distribución para una empresa paraguaya ensambladora de motocicletas, Asunción: Universidad Católica "Nuestra Señora de la Asunción", 2013.