An Immunological Algorithm for Combinatorial Optimization: the Fuel Distribution Problem as Case Study

Artificial Immune systems (AIS) represent a branch of Computational Intelligence (CI), and have been successfully applied to a wide variety of application areas, such as for instance in combinatorial and global optimization [1-3], as well as in systems and synthetic biology [4-6]. AIS are bioinspired algorithms that take their inspiration from the natural immune system and consist of a complex network of interactions among multiple types of agents, whose aim is to detect and contrast all together the antigens, i.e. the foreign organisms that can be cause of pathologies and disease, in order to organism. All AIS algorithms that mimic the clonal selection principle form a special class called Clonal Selection Algorithms (CSA), and today represent an effective mechanism for searching and optimization [7,8]. The core operators of CSA are (i) the cloning, which triggers the growth of a new population of high affinity value, (ii) the hypermutation, which is a search procedure that leads to a faster maturation during the learning phase, and (iii) the aging, whose main purpose is to introduce diversity into the population in order to escape from local optima during the evolutionary search process.


Introduction
Artificial Immune systems (AIS) represent a branch of Computational Intelligence (CI), and have been successfully applied to a wide variety of application areas, such as for instance in combinatorial and global optimization [1][2][3], as well as in systems and synthetic biology [4][5][6]. AIS are bio-inspired algorithms that take their inspiration from the natural immune system and consist of a complex network of interactions among multiple types of agents, whose aim is to detect and contrast all together the antigens, i.e. the foreign organisms that can be cause of pathologies and disease, in order to organism. All AIS algorithms that mimic the clonal selection principle form a special class called Clonal Selection Algorithms (CSA), and today represent an effective mechanism for searching and optimization [7,8]. The core operators of CSA are (i) the cloning, which triggers the growth of a new population of high affinity value, (ii) the hypermutation, which is a search procedure that leads to a faster maturation during the learning phase, and (iii) the aging, whose main purpose is to introduce diversity into the population in order to escape from local optima during the evolutionary search process.
In this research work, we present an immunological algorithm based on the clonal selection prin-ciple designed to tackle a new combinatorial optimization problem that is the Fuel Distribution problem (FDP). FDP is a classical routing problem [9] that a generic fuel transport company faces every-day. Such problem is similar but not identical to the well-known Capacitated Vehicle Routing problem (CVRP) [9], and can be seen also as a variant of the classical Multiple Container Packing problem [10][11][12][13]. The differences between FDP and CVRP are simple but very crucial as they affect both the optimization strategy to use, and the design of the instances to tackle. Given a fixed number of vehicles, each having limited capacity of fuel transportable, the main goal of FDP is to satisfy with the minimal cost all requests made by customers (i.e. fuel stations) that need a fixed quantity of fuel.
It is well known in literature that pure EAs are unsuitable and inefficient when the problem to tackle shows a complex search spaces. However, their hybridization with other techniques might be greatly help for the algorithm. All those EAs that incorporate a deterministic approach, or a local search process in order to refine the solutions, improving then the fitness function, are called Hybrid Algorithms, or better Memetic Algorithms (MAs). The strength point of MAs is the trade-off between the abilities of exploitation by the deterministic or the Randomized Local Search used, and of exploration by the EAs [4,[14][15][16][17][18].
In this research work a clonal selection algorithm is presented for solving the FDP, which incor-porates a deterministic approach for the assignment scheme based on the (DFS) algorithm, and a local search operator based on the exploration of the neighborhood. In order to evaluate the performances of the hybrid clonal selection algorithm, hereafter called H C SA, a real data instance with 82 vertices has been used. Furthermore, in order to understand the robustness and effectiveness of the proposed algorithm we extend the experiments also on 25 different instances, taken from DIMACS graph coloring benchmark, being one of the most popular and used in literature.
The paper is organized as follow: in Sect. §2 we give a general description of the tackled problem with its features and its constraints, presenting in subsection §2.2 a formal definition from a mathematical perspective. In subsection §2.1 we describe why and what are the core differences between FDP and CVRP. In Sect. §3 we present the developed H C SA with its details and main features, dedicating entirely also two sections to the description of the local search designed and the deterministic approach developed for assigning one vertex to one vehicle (Sect. §4 and Sect. §5, respectively). In Sect. §6 we present the study conducted to understand the effects of each parameter with respect to all others, and how it effects on the output, using the Morris method for sensitivity analysis. The experimental results performed and comparisons made are showed and presented in Sect. §7. Finally, Sect.
§8 contain the concluding remarks and the future research perspectives.

Fuel Distribution Problem
The Fuel Distribution Problem (FDP) is a real-world task that any petrol company must face every day in its business activity. FDP is a classical routing problem whose aim is restock along a geographical map a set of customers (usually petrol stations, but not only), minimizing its overall costs. Basically, FDP asks to satisfy all received requests by customers with the minimal cost for the transport company using a fixed number of vehicles, each one having a limited quantity of fuel transportable, called capacity. How to schedule the fuel distribution depends on several questions, as for instance the amount of fuel required in the overall; number of vehicles available to satisfy the demands; relationships among the demands; capacity of each vehicle; road traffic; weather conditions, because the fuel is very sensitive to temperatures (a small amount of fuel evaporates with high temperatures); and many other conditions that affect both road network and amount of the transportable fuel.
Solving such problem means to optimize several objectives subject to some constraints, such as for instance maximize the number of satisfied customers, supplying the exact demand amount by each customer; maximize the number of customers served; minimize the distance from one customer to another; minimize the evaporation along the roads; minimize the costs in each path; and so on. It is important to point out that any feasible route is built taking into account some restrictions, such as (i) the fixed number of vehicles; (ii) any vehicle is able to transport only a limited fuel amount; and (iii) the sum of the amount of a subset of customers satisfied by one vehicle must not exceed the capacity of the vehicle itself.
FDP is then very similar to the Vehicle Routing Problem (VRP), and primarily with its variant better known as Capacitated Vehicle Routing Problem (CVRP) [9]. VRP represents the class of problems in which a set of routes for a fleet of vehicles must be determined for a number of geographically dispersed cities or customers. The objective of the VRP is to deliver a set of customers with known demands on minimum-cost vehicle routes originating and terminating at a depot. CVRP, instead, is the more studied member of the family, where the uniform capacity restrictions for the vehicles are imposed. Such class of problems is closely related to two difficult combinatorial problems: (1) checking whether there exists a feasible solution is an instance of the BPP; (2) setting the vehicles capacity as infinity gets an instance of the Multiple Traveling Salesman problem [19]. Thus, due to the interplay between these two NP-complete problems, and also that the decision version of BPP is conceptually equivalent to a VRP model, the VRP instances can be extremely difficult to be solved in practice [9,19]. Although FDP and CVRP may be seen as the same problem, is important to emphasize that they are not identical due to simple differences that affect crucially the optimization strategy to use.

Differences between FDP and CVRP
In this section we present some crucial differences between FDP and CVRP highlighting and explaining why they are similar but not identical. These differences are crucial, as they affect both the optimization strategy to use for solving the problem, and the design of the instances of the problem. All definitions described below were taken by [9].
(1) In CVRP the given graph ; ) taken in input must be strongly connected and it is generally assumed to be complete, where the cost travel cij to go from vertex i to vertex j, with ( ) , i j E ∈ , is defined as the Euclidean distance between the two nodes. In this way the costs matrix obtained is symmetric and satisfies the "triangle inequality": Satisfying this triangle inequality leads the search process to choose direct links rather than other paths.
In FDP instead is not known a priori what the topology of the given graph is: in general the given graph is a map (e.g. a geographical map) and therefore a sparse graph. It may be a planar graph for small maps, whilst it becomes a more dense graph for bigger ones. Anyway, it never will be represented by a complete graph. Furthermore, since in each edge is associated more than one kind of cost to satisfy (e.g. traffic, distances, weather, travel time, size of the roads), the outcome is the one to not guarantee the property of triangle inequality ( Figure 1). (2) Because in FDP the instance is represented via a sparse graph, it becomes hard or almost impos-sible satisfy one constraint, which instead is one of important conditions for finding feasible solutions in CVRP: each customer vertex is visited by exactly one circuit. Of course, this condition is not true in FDP. Albeit each customer must be satisfied by only one vehicle, more vehicles may instead visit it. This happen because a vertex may be a link point among two or more customers, and thus a crossing among different paths. Besides, each vehicle can cross (then visit) more nodes in its road, not only because the instance isn't a complete graph -becoming almost mandatory to cross more times a given country (example figure. 1) -but also because of the several costs that affecting on the edges that might strongly affect the choice of the path to cross.
(3) Whilst in CVRP are taken into account only vehicles with same capacity, in FDP is not said. Indeed, it is possible to have vehicles with same capacity as well as with different capacities. Of course, this strongly depends by the organizational structure of a given fuel transport company.
These differences, although simple, become crucial for the design and develop of whatever algorithm, as they require different approaches and algorithmic strategies than instead the applied ones in CVRP.

Mathematical formalization as graph theoretic model
Fuel distribution problem can be formulated as follow: let be ( ) an undirected graph, where each node v V ∈ represents a customer (i.e. fuel station), and the set E represents the road network, i.e. each edge ∈ E is the link between two customers' u and v. For the sake of simplification, we assume that In each vertex v is assigned a weight ( ) ( ) : c e c E > →� that represents the costs on the road segment e. In this work we assume that ( ) ) ( c e for all e E ∈ includes all weights that affect the road network (see description above), i.e. the distances among stations (in km for instance); fuel's evaporation degree; temperature; traveling time for each road (hours or minutes), and many others.
Let R = {1, . . . , h} a set of vehicles, such that R V < , and let →� , the limited amount of fuel transportable assigned to the vehicle r ∈ R, called capacity of r. Without loss of generality, we assume that the smallest demand is less than or equal to the smallest available capacity, as well as the largest demand is less than or equal to the largest capacity assigned and besides, that the sum of all weights of vertices in V is greater than the largest capacity assigned to a vehicle in R The aim of FDP is to create h routes in such way to minimize the costs on each route, and maximize the satisfiability of the received demands. For simplicity then we can say that FDP asks basically to assign n vertices to h vehicles, such that (1) the overall number of the assigned vertices to vehicles must be as large as possible (|V | is the optimum); and (2) the sum of the weights of the vertices assigned to each vehicle must not exceed the capacity of vehicle itself. It is worth noting as the constraint (1) means maximizing as much as possible the number of customer demands. Note that this last kind of definition is equivalent to partitioning V in h subsets, such that their union is all V, whilst their intersection is the empty set (θ ).
Although the problem may be seen as a multi-objective problem, we have used instead the following single-objective function in order to evaluate each candidate solution Where V r is the set of all vertices assigned to the vehicle r; β is a penalty value, which assures us to give priority to those solutions able to satisfy all demands of customers; and C tot is the total cost produced by the candidate solution x → , and it is given by Where E r is the subset of all edges visited by the vehicle r.
Furthermore, since the candidate solution x → represents a permutation of vertices from which is determined the visiting order in the graph, it follows that

HCSA -A Hybrid Clonal Selection Algorithm
Hybrid immunological algorithms are now considered a wellestablished technique, as these have been included/studies in [15]. The proposed hybrid immunological algorithm is based on clonal selection principle, whose main features are the operators of cloning, hyper mutation, and aging. HCSA is a modified and adapted version of the well-known opt-IMMALG, already proposed in [2,7,20].
As a typical population-based algorithm, HCSA works with a population of B cells, which has the purpose to defend the organism against foreign, better known as antigens (Ag). From an algorithmic perspective, an Ag is the problem to tackle, whilst the B cell represents a solution for the Ag. In particular, in our study, the Ag is an undirected graph, whose vertices are the customers 1 and edges are 1 i.e. fuel stations that have a contract with the distribution company, included the ones that haven't made demands the road connections between customers (see Sect. §2.2). Each B cell, instead, represents a permutation of vertices (i.e. a string of integers) that defines the visiting order of the customers for the assignment to a proper vehicle (if it exists) ( Table 1).
In table 1 is presented the pseudocode of HCSA, where by P (t) is indicated a population of d individuals of length l = |V |. The initial population is randomly generated by a uniform distribution, which represents a subset of the search space, i.e. the starting points for the searching in the space of solutions (line 1, pseudocode). How to generate the initial population is a crucial task for EAs, as it might influence the overall performances. In traditional EAs, the initial population is generated using a random numbers distribution or chaotic sequences [21]. An interesting research work, which studied the properties of different point generators for a genetic algorithm, was presented in [22].
Miming the clonal selection principle, HCSA incorporates the static cloning operator, that clones all B cells dup times, producing an intermediate population P (clo) (line 5, pseudocode). In this work we haven't used a proportionally cloning, as it shows premature convergences. This is due to a strong pro-liferation of the best solutions found so far, which decrease the diversity in the population, and therefore not help H C SA to escape from local optima. Afterwards, each cloned B cell is undergo to the hyper-mutation operator, which mutates the cloned M times, without an explicit usage of mutation probability (line 6, pseudocode).
of mutations is determined in an inversely proportional way to the fitness value, although there exists several approaches as proposed in [23]. The mutation rate is determined by where f is the fitness function normalized in the range [0, 1], and ρ is a parameter. The Figure 2 shows the curves produced by the equation 3 at different ρ values (left plots), and the number of mutations obtained for different problem dimensions, and ρ values (right plot). For each B cell receptor, the hypermutation operator chooses randomly times two vertices u and v in x → , and then it swaps them having as effect to change the visiting order of the graph. At the end of the hypermutation process, we have a new population that is denoted by P (hyp) . In order to normalize the fitness function in the range [0, 1], as proposed in [3,7], HCSA uses the best current fitness value decreased of an user-defined threshold Θ ; this is due because is not known a priori any kind of information about global optima. In this way, HCSA doesn't need to have any kind of information concerning the problem; HCSA is a real blind box. In this work we have set Θ = 50%. One crucial question that might affect the performances of any immunological algorithm is what age assigns to each hypermutated clone. In order to give an enough evolutionary time for the maturation for each B cell, HCSA uses an equal opportunity scheme: when a hypermutated B cell improves the value of the fitness (called constructive mutations), then it will be considered to have age equal to 0; otherwise, it will maintain the same age of its parent. By this scheme we want to give an equal opportunity to each new B cell to effectively explore the search space. To refine the quality of the solutions, HCSA incorporates also an heuristic local search based on the exploration of the neighbourhood, where the neighbours are generated through the swapping of the vertices. This operator is described in section 4, and it is applied to the best receptor of P (hyp) , producing a new population, called P (LS) (line 8, pseudocode).
After the perturbation operators, all the old B cells inside the populations P (t) , P (hyp) , and P (LS) , are eliminated by the static aging operator (line 10, pseudocode). The parameter τ B in input indicates the maximum number of generations, that allows to each B cell to remain into the corresponding population: when a B cell is τ B + 1 old it is erased from the own population, independently from its fitness value. This means that each B cell is allowed to remain into the population for a fixed number of generations. An exception is made only for the B cell with the best fitness value (elitist static aging operator). The aim of the aging operator is to produce a high diversity into the current population, and avoid premature convergences. Therefore, the aging operator plays a central role on the performances of the proposed algorithm (and AIS in general): too much diversity inside the population could produce poor solutions, as well as too small ones.     new B cells will be randomly created. (Figure 3) Finally, Evaluate _Fitness(P ( * ) ) (lines 2, 7 and 9, pseudocode) computes the fitness function value of each B cell x P → ∈ using Eq. 1, whilst Termination Condition() (line 4, pseudocode) is a Boolean function, which returns true if the maximum number of generations, or the maximum number of fitness function evaluations allowed, is reached; false otherwise. It's important to note that in order to minimize the costs of the routes, HCSA includes the Dijkstra algorithm [24], which computes for each vehicle the shortest path between two unconnected customers u and v, that is when ( ) , u v E ∉

Local Search
As cited in the previous section, afterwards the hypermutation operator, HCSA incorporates a heuristic local search in order to refine and improve the quality of the solutions. The used approach for the design of the local search was taken from [25,1], and it relies on the definition of neighbourhood, where neighbours are generated through the swapping of the vertices.
Starting from the idea that the visiting order of a graph shapes the solutions found, we define the neighbour of x → , all B cells y → that can be obtained from x → by swapping two of its elements. Because swapping all pairs of vertices is time consuming, we have used a reduced neighbourhood by a radius R LS , as proposed in [1,2]: in each B cell, all vertices were swapped only with their R LS nearest neighbours, to the left and to the right. Moreover, take into account the large size of the neighbourhood, we have applied the local search procedure only on the best hypermutated B cell (i.e. the best of P (hyp) ). When a single swap between two genes reduces the fitness function value (constructive mutation), then the new mutated B cell is added into the new population P (LS) ; otherwise it is not taken into account, and hence erased. The process continues until the whole neighbourhood with radius R LS is explored. To avoid the problem to study what is the best tuning for the R LS radius, HCSA randomly assigns a value in the range [1, (|V | − 1)], using a uniform distribution [25]. In this way it guarantees to swap at least two vertices. Figure 3 shows the curves of the evolution of the best fitness using the local search procedure, or no. The experimental was made fixing the minimal values for the parameters: d = 100, dup = 1, τ B = 5, ρ = 5.5 and M axGen = 1000. From the figure is possible to see how the use of the local search helps the convergence process towards better solutions.

Heuristics for the Assignment Scheme
How to assign a vehicle to one customer or vice versa is a central point in the design of the algorithm. Because each vehicle has a maximum capacity of transportable fuel, then choosing one vertex, rather than another, can determine the satisfiability of all demands (the goal), or only some of them. To this purpose, in this work all B cells represent a vertices permutation, which determines the visit order. In details, the used heuristic scheme works as follows. , is possible to distinguish the following cases:

if exists a vertex v V
∈ adjacent to x j , and a vehicle r R ∈ assigned to v, such that where V r is the subset of the vertices already assigned to the vehicle r, then r is assigned to x j . If there exist two or more vertices adjacent to x j , with assigned different vehicles suitable to satisfy x j , then one with higher available capacity is assigned; 2. for all vertices v V ∈ adjacent to x j , or not exist any r R ∈ assigned to v, or if there is, it is not able to satisfy x j . Thus, if there are one or more free vehicles, i.e. not still assigned, then one of these is randomly chosen, and assigned to x j otherwise a deep search is made into the neighbourhood of x j . If after the search, at least one vehicle was found, this is assigned to x j , otherwise the vertex will be labelled "not satisfied". In this work, as search model, was used the classical "depth first search" algorithm (DFS) [24], but reduced of a radius R (DFS) : is fixed a limit R (DFS) < |V | to the depth of the search into the neighbourhood. The purpose of this reduced DFS is to find a suitable vehicle not too far from x j , in such way to have homogeneous groups. If, by the reduced dfs, we found two or more suitable vehicles, then the nearest one is assigned to the vertex x j . For the experiments described in section 7, R (DFS) was fixed to 15% of V. Such value was chosen to avoid solutions where one vehicle must satisfy two stations too far away, for example geographically opposite.
3. not exists any vehicle able to satisfy the given vertex: i.e. for all where V r is the subset of the vertices assigned to the vehicle r. In this case the vertex will be labelled "not satisfied". Of course this occurs when the two previous steps fail.
After any assignment, the capacity of the chosen vehicle r is decreased: b curr (r) = b prev (r) − q(x j ). With regard to the approach 2, if there exist more vehicles still unassigned, then one of these will be randomly chosen for the assignment, and this scheme is called "random + dfs assignment". Another variant may be taken into account, that we call "dfs 2 assignment": before randomly choosing a free vehicle, this variant checks the entire for a neighbourhood suitable vehicle that could be assigned; the reduced dfs, where the radius R (DFS) is fixed to 5% of V , is used. This new scheme, guarantees us the design of more  Figure 4 shows the comparisons of the curves of the best fitness obtained by the two described approaches. These curves were obtained using the following parameters: d = 100, dup = 2, τB = 15, ρ = 5.5, and M axGen = 1000. This figure shows how the basic idea designed in dfs 2 allows HCSA a better convergence towards solutions with best quality, and lower costs. The inset plot shows the curves of success of the DFS (best run) for the both approaches. Obviously, as we expected, the curves of random + dfs are higher than dfs 2 , because it develops solutions more poor, and then it needs more ( Figure 4) calls to the long DFS, i.e. with radius R (DFS) = 15% of |V |. Therefore, more calls to DFS, generates an higher number of DFS success. However, the overall percentages of the DFS success, averaged with their calls and computed on 10 independent runs, is higher in dfs 2 (73.08%) than in random + dfs (66.28%). In table 2, instead, is showed the comparison between the two cited approaches, varying the parameters, on the graph with 82 vertices, based on real data. This instance was used to evaluate the performances of HCSA not only with respect the fitness value, but also and primarily in the develop of homogeneous groups. In this table is shown the best fitness values, the mean and the standard deviation (Table 2).
(σ). Last column ∆ of the table indicates the differences of the two approaches with respect the best fitness values. In bold face is highlighted the best fitness values for each pairs of parameter. The results were obtained with d = 100, ρ = 2, M axGen = 100, and 10 independent runs. The number of the vehicles h was fixed to 6, using the same capacity for all vehicles. Given the set of weights on the vertices, we can compute as an obvious lower bound (but not optimal) on the vehicle capacity needed to satisfy all required by the stations. In this work, we have fixed as capacity of the vehicles the lower bound λ increased by 0.2%; i.e. as small as possible. Inspecting such table is possible to see how dfs 2 is more suitable in finding better solutions, as well as more homogeneous. The better performances produced by dfs 2 are due because after the first generations the called number of the long DFS, i.e. the ones with radius R (DFS) = 15% of |V|, decreases, having developed current subgroups homogeneous. Therefore, will be easier to find a suitable vehicle into the near-neighbourhood. Thanks to that, dfs 2 is more faster in computational time: the computational times are respectively 16m25919s for random + dfs, and 8m17953s for dfs 2 , obtained using the parameters d = 1000, dup = 5, τ B = 5, ρ = 5.50, and gen = 200.

Sensitivity Analysis
In sensitivity analysis methods, one can observe how a parameter affects the complexity of the instance, while all other parameters are varied simultaneously. These methods consider the interactions between parameters without depending on the stipulation of a nominal point.
The where ∆ is a predetermined multiple of 1/(k − 1). The distribution of elementary effects F j is obtained by randomly sampling h points from Ψ.
For each parameter we evaluate two sensitivity measures, such as (1) µ j an estimate of the mean of the distribution F j , (2) and σ j an estimate of the standard deviation of F j . A high value of µ j indicates a parameter with an important overall influence on the output; whilst a high value of σ j indicates a parameter involved in interaction with other parameters or whose effect is nonlinear. In order to understand what are the parameters that affect the output, we performed the Morris method on a graph with 256 vertices, and 32640 edges (queen16 16.col -see section 9). The Morris method is a sensitivity analysis useful to understand the effects of a parameter with respect to all others, which vary simultaneously. Figure 5 shows the sensitivity analysis carried out for our objective function (equation 1). Inspecting this figure is possible to see how the vertices {36, 115, 139, 152, 172, 234} seem to be the most important, since they affect more on the objective function than the remaining vertices, whereas nodes 118, 101, and 182 are the less influential ones. (Figure 5)

Results
To evaluate the goodness of the performance of HCSA and its search ability into the solutions space, we have used two different evaluation measures: (1) if HCSA is able to obtain good approximate solutions using the capacity of the vehicles as small as possible, and (2) the homogeneity in the assignment of the vehicles to the vertices; i.e., to avoid that a vehicle has to supply two vertices placed in opposite sites from a topological point of view. For the experiments, we have used initially a graph, based on real values, with 82 vertices. Afterward, to extend our test beds we have tested HCSA on several graphs, taken by the Dimacs colouring benchmark [27], being one of the most popular and used in literature. Therefore, in this section we report all results obtained in our experiments, describing also the experimental protocol used for each test. Of course, once experimentally proved that dfs 2 shows better performances than random + dfs with respect to the costs and the homogeneity of the solutions (see table 2), then, all results presented in this section were obtained by the dfs 2 heuristic. In table 3 we report the results obtained by HCSA using for all vehicles either the same quantity of the transportable fuel, and different capacity. The capacity was increased by 0.2% of the lower bound λ -Eq. 5. These experiments were made varying the parameters dup = {1, 5, 10}, and τ B = {5, 20,..., ∞}, and fixing population size d = 100, ρ = 4 for the experiments where all vehicles have the same capacity, and ρ = 5.5 for all experiments with different capacity values. Moreover, the maximum number of generations was fixed to 100, and for each test was made 10 independent runs. If we give a look to the results obtained using the same capacity for all vehicles, one can see that, although the best solution is obtained with high values of dup (dup = 10), in the overall the better performances are instead obtained using smaller values; this is proved by the mean values. If we use different capacities for the vehicles, then dup = 5 seems instead to be the adequate setting. To simulate a real world application, we have considered the graph as road network, where each weight has been randomly generated in the range [200, 10000] (as litres) for the vertices, while for the edges in the range [5,180] (as minutes). Moreover, we have used a small (Table 3) number of providers in order to better simulate a real application. Since in the real world is unlikely that all customers of a distribution company make demand at the same time, i.e. some node can have weight null, the random generator assigns each weight on vertices with a probability P = 50%.
Understanding the real capabilities by HCSA from an exploration and exploitation perspective, we have compared HCSA with a classical Genetic Algorithm (GA) and a deterministic algorithm based on locally optima choice strategy. For this deterministic algorithm we present three different versions: (1) starting from the vertex V1 the naive method proceeds sequentially (V 1 , V 2 , ..., V n ). We labelled this version as naive; (2) starting from a random vertex Vk , this method proceeds as follow (V k , V k+1 ,..., V n , V 1 , ..., V k−1 ). We call this second version as DBO; (3) the last method performs the optimal locally selection based on a permutation of the vertices (DBP). Table 4  A high value of _indicates a parameter with an important overall influence on the output. A high value of _ indicates a parameter involved in interaction with other parameters or whose effect is nonlinear.  [1], being one of the most popularin literature. For each instance is showed the number of items satisfied (Γ), and relative best cost found. The experiments have been performed for 30 independent runs. We point out that if one of the algorithms is not able to satisfy all requested of the items, the relative costs have been not included in the table (_).
stop criterion a maximum number of fitness function evaluations T max = 5 × 10 4 for all graphs with |V | < 100, and T max = 5 × 10 5 otherwise. Besides, 30 independent runs have been performed (Table 4) for each instance. H C SA is able to satisfy all demands received, with respect deterministic algorithms. HCSA and GA have been able to satisfy all request received on different dimensions of the problem (from 36 to 256 vertices). However, comparing HCSA and GA is possible to see how the proposed algorithm is able to find better costs in all tested instances, which means that HCSA is able to produce more compact groups from a topological point of view. In the table 4, if one of the algorithms has not been able to satisfy all requested, the relative costs ( f ( x → )) have been not included in the table and labeled as −, because the fitness value produced is high due to the penalty factor β (Eq. 1).

Conclusions and Future Work
In this research work we present an hybrid clonal selection algorithm (HCSA) for the Fuel Distribution Problem, one of the classical combinatorial optimization problem that many transportation companies must face in their everyday operations, that is to restock all their customers along a geographical map, minimizing their overall costs. Such problem is very similar but not identical to the Capacitated Vehicle Routing Problem -CVRP, and can be also seen as a variant of the classical Multiple Container Packing Problem -MCPP. We present in a proper section (Sect. §4.1) the simple differences between FDP and CVRP, which significantly effect on the optimization strategy to use, and on the design of the instances to tackle.
HCSA is based on three main operators: (i) the cloning that triggers the growth of a new population of high-value B cells; (ii) the hypermutation that can be seen as a search procedure that leads to a fast maturation; and (iii) the aging that causes a turn-over in the populations with the aim to generate diversity and avoid to get trapped into local optima. HCSA incorporates a local search heuristic in order to refine the solutions found so far, and it is based on the exploration of the neighbourhood, where the neighbours are generated through the swapping of the vertices. This heuristic is founded on the idea that the visiting order of vertices significantly affects the satisfiability or less of all demands made by customers. For proving the usefulness and improvements produced by designed local search, a study on the impact of such local search has been conducted. Such study has confirmed us how such heuristic helps HCSA in refine the quality of the solutions, and giving a best convergence, towards best solutions. In this research work, we propose also two variants for the assignment the customers to the most suitable vehicle, both based on the DFS algorithm. The first variant is a combination between random assignment, and DFS algorithm; whilst the second, before to choose a free vehicle for the assignment, the scheme performs a search in the neighbourhood of radius R (DFS) . In order to limit the length of the search tree, we applied a reduced DFS based on a radius R (DFS) ; this reduction is to avoid solutions where one vehicle could be assigned to two vertices too far away. From the conducted experiments the reduced DFS seems develop good solutions, in term of quality, and homogeneity of assignments. It seems to be very promising for the future direction of our research on this problem.
To evaluate the performances of HCSA, we have used a graph with 82 vertices based on real data. However, to extend our experiments, and comparisons, we have also used some graphs with different topologies that were taken by the Dimacs graph colouring benchmark, being one of the most popular and used in literature. For this last class of experiments, we have randomly generated the weights on the vertices, and on the edges. All results were obtained using as capacity for the vehicles the lower bound λ (Eq. 5) increased by 0.2%, i.e. as small as possible, to really understand the search ability of HCSA into the solutions space. Furthermore, HCSA has been compared with a Genetic Algorithm, and three different versions of a deterministic algorithm. Inspecting the results is possible to clearly see how HCSA outperforms the compared algorithms in all instances. HCSA is able to satisfy always all demands, as opposed to the compared algorithms, which instead fail in several instances.
Several new research directions need to be investigated: (1) study the best tuning of the radius R (DFS) for the reduced DFS; (2) design a better refinement operator, such to improve the convergence speed of HCSA, and its computational efforts; (3) and finally, perform and test HCSA taking into account a dynamical environment where the weights on the vertices and/or on the edges may change during the time, as well as vertices and /or edges turn on or turn off in the time.