Path relinking based intensification strategies for a simulation-optimization scheduling problem arising in hydroinformatics

Abstract

Water contamination events in a urban hydraulic network can be handled by operating two types of network devices: valves, to divert flow, and hydrants, to dispel flow. No remote control can activate these devices, which must be operated on site, therefore teams of technicians have to move from device to device on the city streets. The response to the contamination corresponds to a feasible schedule of the devices’ activation times, which can be modelled as a multiple Travelling Salesman Problem; the quality of the schedule is the volume of contaminated water consumed by the users until water quality returns to safety levels. This quantity is highly sensitive to the schedule and is computed by way of a computational demanding hydraulic simulation, which takes about 5 seconds for the Ferrara’s hydraulic network serving around 130’000 citizens. Previous studies proved that minimizing the makespan, as it would be intuitive in case of emergency, yields worse solutions than random approaches. Genetic Algorithms were proposed to optimize the schedules, however they converge to local optima, depending on the initial population. We propose to apply Path Relinking to explore the space between pairs of local optima; the experiments showed that this technique improves the local best. This work is a follow up of what presented at AIxIA 2015, enpowered with a light introduction to Path Relinking and its applications. Moreover, the experimental campaign has been extended on a larger set of scenarios, i.e. 42 contamination scenarios, to gain more insight into the performance of the proposed methodology.

Keywords

Scheduling path relinking simulation-optimization hydroinformatics

1 Introduction

Hydroinformatics is an interdisciplinary research field arising at the junction of Hydraulic Engineering and Computer Science, in which complex decision problems related to water management applications are modelled and solved by quantitative solution tools developed within well assessed computer science methodological paradigms, such as Constrained Programming on Finite Domain and Mathematical Programming Optimization. Several such examples can be found in the literature which exploit the network based problem structure, taking advantage of solution methodologies already developed for transportation and communication networks, as earlier pointed out by Simonis in a seminal work [41]. Among the most recent contributions, let us mention: the design of the expansion of a Water Distribution System (WDS) combining global search techniques with local search [5], the optimal location on the WDS of water quality sensors in order to early detect water contamination exploiting integer programming based location models [30] and objective function sub-modularity [25], the optimal scheduling of devices controlling field irrigation [22] to meet farmers irrigation time demands, the optimal location of isolation valves on the pipes of a WDS to minimize service disruption in case of maintenance operations [8] and [11], and the scheduling of devices activation as a countermeasure to contamination events [12], which is the problem we deal with in this paper.

In all cases, the feasible solutions have to meet complex technological requirements which are modelled by the constraints of the optimization problem, while the objective function describes how the hydraulic system reacts to certain values of the parameters, which are the model variables. Quite often, the system reaction can not be encoded by analytical closed formulas but it is the result of a computational demanding simulation process, which poses a challenge to the development of solution methodologies able to scale efficiently and tackle real life instances.

In this paper we deal with the last mentioned problem, namely computing the optimal activation time of a set of devices located on the WDS. In case of a contamination event, the devices activation times alter water flow, influence how contaminant spreads in the network, and determine at which concentration contaminant reaches demand nodes where drinking water is consumed by system users. An optimal schedule is a set of activation times which minimizes the volume of consumed contaminated water. A feasible schedule is a set of activation times according to which the teams of technicians, who manually activate the devices on site, can reach the selected device on time, travelling on the road network when moving from a device to the next. Common practice at most water distribution agencies is to compute such schedules by hand, driven by common sense principles, such as “the sooner the better”. The alarm is raised by the activation of one or more sensors. Since there are only a few, it makes sense to enumerate a representative set of contamination scenarios, differing on i) the activated sensors, ii) the location and iii) the time of the pollutant injection. For each scenario the activation plan is computed and stored. For unexpected scenarios there is no time for computation and either common reactions or the precomputed plan for the closest scenario are the only available. Notice that for experimental needs we test only 42 scenarios, however the methodology can be applied to any expected scenario. Recent approaches exploited genetic algorithms (GA) to compute better schedules automatically [2 , 32]. In this paper, we extend the study presented at AIxIA 2015 [31] devoted to investigate how hybrid solution approaches for such complicated real problems can largely benefit from the integration of different search paradigms. We also believe that path relinking, thanks to its applicability to basically any combinatorial optimization problem, could be as well integrated within most CP engines, acting as an intensification module to be activated on any set of good quality (partial) solutions visited so far in the search tree.

In the rest of the paper, first we introduce the problem and recall the solution strategy based on Mixed Integer Linear Programming (MILP) hybridized genetic algorithms developed so far, pointing out at some deficiencies. Then, we describe a neighbourhood based search strategy, so called Path Relinking (PR) originally proposed by Glover [18], which intensifies the search within a section of the feasible region to which a set of high quality solutions belong. We present how to use PR in our problem, compare two different neighbourhood structures to build the search path connecting two solutions, and assess the effectiveness of the approach by experimental results computed on real data for the WDS of a medium size city in Italy, showing how by enhancing our GA with a post optimization PR phase can improve the approach robustness and partially mend the present flaws.

2 Problem description

WDSs are essential components of our daily life as they bring clean, safe drinking water to customers every day. At the same time, WDSs are among the most vulnerable infrastructures, highly exposed to the risk of contamination by chemical and biological agents, either accidental or intentional. A WDS is a complex arrangement of interconnected pipes, pumps, tanks, hydrants and valves, whose large planimetric extent (a small city network may reach 200 km and a thousand of pipes and nodes) and sparse topology make full surveillance impossible. Therefore monitoring is the only viable alternative. In practice, a set of water quality sensors is located at different sites on the WDS to test water safety in real time, looking for the presence of potential contaminants [30]. Their location is strategically determined so that acontamination event is detected as soon as possible, based on a set of contamination scenarios.

Contaminant quickly spreads through the network and population alerting strategies may not entirely ward off users’ water consumption. When the network is fully districted, the sector where the alarm is raised can be seamlessly disconnected from the rest of the network, but this is rarely the case. In general, despite of the hazard, water supply can not be completely cut off. The shut down of the entire system would disrupt those security related functions that rely on continuous water supply, such as fire police service or water based cooling systems atlarge, production intensive, industrial facilities. Therefore, beside population warning procedures, countermeasures devoted to divert the contaminant flow away from high demand concentration sectors must be set up, aiming at mitigating population harm.

An effective way of altering water flow is by activating some of the devices which are part of the system, namely by closing isolation valves and opening hydrants, in order to achieve contaminant isolation, containment, and flushing. In particular, opening hydrants can expel contaminated water, while contaminated pipes can be isolated by closing their isolation valves. Due to the highly non linear functional dependencies that link water flow and the time at which a given device is operated, the global effect of a schedule, i.e., the vector of activation times for the selected devices, can not be decomposed into the sum of the effects of the activation of each individual device, nor the effect of a local change in the schedule can be anticipated. On the contrary, the only way to evaluate the volume of consumed contaminated water due to a schedule is by a computationally intensive simulation. In other words, we are optimizing a black box function and solving a so called simulation-optimization problem [3]. The chosen simulation package is EPANET [39], a discrete event-based simulator which represents the state of the art in Hydraulic Engineering. EPANET is given a schedule, the description of the hydraulic network, the expected water demand over time, and a contamination scenario, and it yields the volume of contaminated consumed water (we speak of contamination when concentration level is above the danger threshold of 0.3 mg/ml that causes human death if ingested). Simulation is computationally intensive and it is the bottleneck of any solutionapproach.

In our case study, each simulation call takes on average 5 seconds, which poses a limit on the number of function evaluation calls, despite of the fact that the problem is solved offline; this fact influences the choice of the search strategy, as discussed in [13]. Moreover, there is no a priori information on what a good schedule should look like, and common sense inspired criteria such as “the sooner the better” lead to low quality solutions. In conclusion, there is no way to distinguish between a good and a bad schedule without simulation.

So far we have discussed the features of the objective function of our simulation optimization problem. Regarding solution feasibility, a schedule $t^{F}$ is feasible provided that there is an assignment of the n devices to the m teams available and, for each team, its devices can be sequenced so that if device j follows device i the difference between the respective activation times in the schedule is equal to τ_ij, i.e., the travelling time on the street network from the location of device i to the location of device j. All teams gather at the mobilization point at a given time after the alarm is raised (according to the protocol) and conventionally the departure time is set to 0. This maps the feasible region of our problem into the one of a well known combinatorial optimization problem, the open m-Travelling Salesman Problem [4] (mTSP). This provides a graph representation of our problem where the mobilization point is the depot, each device is a client node of the graph and each team visits the assigned devices travelling along a route starting from the depot and ending at its last visited client.

The resulting mathematical model for this simulation-optimization problem is given as in [13, 32]. The input parameters are:

a set of devices 1..n.

a matrix (n + 1) × (n + 1) (including the depot); τ_ij represents the time that a team takes to move from the location of device i to that of device j; a given constant speed of the teams is considered on the street layer of the city.

The unknowns are:

a matrix (n + 1) × (n + 1) of 0-1 variables. x_ij = 1 iff j is activated right after i by the same team; i is activated first by its team iff x_di = 1; x_ii = 0 ∀i (no self loop arcs).

a vector of n + 1 activation times; $t_{i}^{F}$ is the time at which device i is activated, and $t_{d}^{F}$ is the departure time from the depot d.

The constraints are: $t_{d}^{F} = 0$ (1a) $\sum_{i \in D} x_{di} = m$ (1b) $\sum_{j \in D \cup d} x_{ij} = 1, \forall i \in D$ (1c) $\sum_{j \in D \cup d} x_{ij} = \sum_{h \in D \cup d} x_{hi}, \forall i \in D$ (1d) $t_{i}^{F} \leq M + x_{di} (τ_{di} - M), \forall i \in D$ (1e) $\begin{matrix} t_{i}^{F} + τ_{ij} \leq t_{j}^{F} + (1 - x_{ij}) M + \\ + x_{ji} (τ_{ij} + τ_{ji} - M), & \forall i, j \in D \end{matrix}$ (1f)

The aim is to minimize $v (t^{F})$ subject to (1a-1f). $v (t^{F})$ is the volume of contaminated water consumed by the users, and it is computed by the hydraulic simulator. Teams leave the depot at time 0 (Equation 1a)). Each team departs from the depot (Equation (1b)). Each node except d is visited exactly once (Equation (1c)). For each node i, the total number of teams arriving to i is equal to the number of teams leaving i. Equation (1d) are the so called flow balance constraints. Constraint (1e) is the linearisation of the implication $x_{d i} = 1 \Rightarrow t_{i}^{f} \leq τ_{d i},$ where M is a sufficiently large positive number. Constraint (1f) links the activation times $t^{F}$ to the ordering between devices given by matrix X.

This special structure leads us to set up a genetic algorithm hybridized with a Mixed Integer Linear Programming solver which encapsulates the feasibility constraints within the cross over operators, as described in detail in [13]. In that paper, several computational experiments were carried out to calibrate population size and number of generations for a single GA run. Moreover, we verified the poor quality of the solutions obtained according to heuristic criteria, such as minimizing the makespan or minimizing the sum of the activation times. Besides, neighbourhood based local search were tested and proved to be not competitive given the limited number of solution evaluations allowed, since the search trajectory remains confined not far from the starting point. On the contrary, the literature confirms that in such cases population based heuristics, which carry on a broader search and are able to explore a wider part of the feasible region, are able to provide better results [3].

Although we could improve by far and large the best solutions available in the literature for our case study, that solution approach has a typical GA flaw, that is, it converges to a local optimum which depends on the starting population. However, the differences among different local optima, obtained by repeated runs starting from different initial populations, vary depending on the contamination scenario. Since the number of function evaluations is limited, in this study we address the question concerning what is the best way of exploiting the available computational resources, and if there is a way to take advantage of the knowledge of a set of different, high quality solutions. Previous studies proved that for these problems population-based evolutionary approaches are most effective [3, 13]; for this reason, in this study we investigate the effectiveness of a well-known evolutionary intensification strategy in the literature, namely the Path Relinking (PR). A comparison among PR and other intensificationstrategies, e.g., by GAs, is under study.

In the next section we first provide a light introduction to Path Relinking and then we describe how it can provide a pattern search that fulfills these expectations.

3 Intensification by path relinking

As mentioned, in this application GAs often yield high quality solutions that depend on the initial population: this is due to the existence of several local optima. These different solutions identify a promising subregion of the feasible space, which is worth further inspection according to some exploration strategy. This paper is devoted to investigate the potential of Path Relinking as an intensification strategy in this specific setting. We believe that the Path Relinking search paradigm fits well into this framework where few solution evaluations are allowed but at the same time there is a set of local optima worth to be exploited as a whole. We first provide a light introduction to Path Relinking fundamentals and recall its applications to routing problems. Then we describe two alternative implementations of the PR paradigm for our problem, based on two different neighborhoods.

3.1 Path relinking fundamentals

PR was first introduced as a complementary component of a new, multi-thread, evolutionary, search strategy called Scatter Search (SS) [18]. SS shares with GAs the idea of generating new solutions by combining current solutions, however the way SS combines them is substantially different from GA since SS does not rely on randomness but it is rather deterministic. Basically, PR concerns the way one or more solutions are generated from other solutions, while SS concerns the policies according to which the entire set of current solutions evolves.

SS roots can be traced back from a seminal paper by Glover [17] describing a way of computing suboptimal solutions to integer linear programming problems by recombining problem constraints via linear combinations. This search pattern evolved over time, adapting its fundamentals to encompass combinatorial optimization problems. In a few words, SS works on a limited set of elite solutions, usually ≤ 20, denoted as the reference set (rs), where diversity among the members is required as well as a good quality of the objective function value. Improvement methods are applied to individual solutions as soon as they are generated, followed by an updating phase of the elite solution set where new solutions are allowed to enter the rs and other solutions are discarded. Then, subsets of solutions (generally two solutions) are selected to be combined according to the PR strategy. Intensification and diversification mechanisms usually enrich the selection step.

It clearly appears that SS shares many components with several other search methods. In particular, the ties of SS with other metaheuristics are analyzed in [6].

PR explores the solutions encountered along a trajectory connecting a pair of given solutions; these are usually referred to as the initial (r) and the guiding or target (g) solution. This search pattern is so general that PR has become a popular and successful intensification strategy on its own, apart from SS, and can be applied whenever a potential reference set is available. Its description requires to formalize a solution s belonging to the feasible set S as a vector of n elements, where element s [i] is the value of the ith attribute of solution s. Now consider a neighborhood function N defined on S, where N (s) is the set of solutions s′ ∈ S such that s′ is sufficiently close to s, meaning that most attributes have the same value in s′ and in s. Generally, N is operationally defined by describing which attributes of s can change in s′ and how they can be changed. Given N, let us introduce a graph G = (S, E) where there is a node for each solution s and there is an edge connecting pairs of solutions that belong to each other’s neighborhood; formally (u, v) ∈ E ↔ u ∈ N (v) and v ∈ N (u). PR considers a special neighborhood structure dependent on g, made of all solutions obtained by changing the value of one attribute of the current solution u into the value such attribute has in the target solution g. In this sense, one often speaks of a restricted neighborhood. In basic PR, the search starts at r, and at each iteration the search moves forward from the current solution u to v^* ∈ N (u), the best solution (according to the objective function) among the adjacent ones in the restricted neighborhood, making the next solution one step closer to g. PR returns the best solution along this path. The number of iterations depends on how much different r and g are, which explains how crucial the issue of diversity in the reference set is. In particular, if r differs from g for a single attribute, the search path is made of a single edge and no new solutions will be explored. Likewise, if both are local optima and they differ for two attributes, then the only solution in the search path does not improve neither on r nor on g. Note that, for a given problem, each alternative solution representation in terms of attributes gives raise to a different neighborhood structure and therefore it provides a different search trajectory to be inspected. We will indeed present two alternative PR implementations, each one associated to a different neighborhood structure.

PR differs from classical Local Search (LS) in many aspects: LS always converges to a local optimum after an unknown number of steps; PR converges to g after a number of steps bounded by n and given by the number of attributes i for which r [i] ≠ g [i], i.e., the distance from r and g. At each iteration LS potentially explores all solutions in N (u) and their number tends to be constant along the search; in PR, N (u) contains at most as many solutions as the number of attributes i such that u [i] ≠ g [i] and this number decreases the closer the search gets to target. In both cases, though, the search moves on the best solution in the neighborhood.

Figure 1 shows a graphical representation of the transformation of r into g. r differs from g by 4 elements; therefore, at the first iteration the neighborhood is made of 4 solutions, these are 3 at the second iteration, and so on; the search takes 4 iterations during which it moves along 3 intermediate solutions, namely r¹, r², and r³, the best of which is returned. Note that the neighborhood of the initial solution is explored more thoroughly than the neighborhood of the final one. This has motivated alternative strategies according to which the search starts from the best solution and moves towards another solution (backward PR), or it searches the path in both directions (back-and-forward PR). None proved to dominate the other strategy.

Another variant is obtained by adding some randomness in the process, so far entirely deterministic: in the same way as Greedy Randomized Adaptive Search [9] extends a greedy solution to a set of solutions by randomly operating locally suboptimal choices, PR can be made exploring several alternative paths linking r to g by choosing at each iteration a suboptimal solution in the restricted neighborhood. By repeating this process a few times, alternative search trajectories can be inspected, further intensifying the search in the area.

Since PR builds a new solution starting from the features of two elite solutions, PR can be also seen as an evolutionary algorithm, in which randomness is substituted by a deterministic search strategy that draws the possible path between two feasible solutions. Reversely, GA can be seen as a randomized blind PR in which one solution along the path is randomly selected and picked, ignoring all the others. This fact confirms the idea that PR should be rather used as a post optimization, intensification procedure once that a promising region has been identified in terms of the solutions in the reference set, for example resulting from GA. However, the two techniques can be hybridized in a completely different way, as in [38]: given two parents solutions, a back-and-forward PR search is carried out and the solution is considered as the result of a crossover operation. For all these reasons PR appears to be a promising intensification strategy when coupled with a GA, as in our case.

The building blocks of a path relinking algorithm are:

the reference set and the criteria for building it;

the selection of initial solution and target solution;

the path between two solutions, i.e., the neighborhood structure.

Several variants and generalizations of PR can be obtained by different implementations of the building blocks, as widely discussed in [19], such as truncating the search on a path to resume it on another (either new or existing) path of the same g – r couple, or different policies for choosing the move in the current neighborhood rather than moving on the best one. All these options involve different trade-offs between solution quality and computational resources. In this work we adopt this last classical strategy since we prefer to spend the limited number of solution evaluations we are allowed to explore the "best" path between different r – g pairs rather than several paths between the same r – g pair.

In our case, the set of best populations given by some GA runs naturally provides the dataset on which the reference set can be built, the target g can be selected as the best solution in rs, and the reference r has to be selected properly through quality and diversity criteria, as Section 3.3 reports.

Note that the quality of the reference set affects the path relinking effectiveness because its search is within the area identified by the reference set. Thus, it can be seen experimentally that the worse is the reference set the lower is the effectiveness of the approach.

3.2 Path relinking in routing problems

From the above discussion, one might deduce that applying PR to a problem for which a neighborhood based heuristic is available should be a straightforward task, since that heuristic could be used to compute a set of reference solutions and PR could make use of the same neighborhood structure N used in the search. This is only apparently true, since there may not be any guarantee that g and r are connected in the graph induced by the restricted neighborhood. Moreover, one may be obliged to arbitrarily establish some correspondence between g and r to be maintained for all the intermediate solutions inspected during the search, in order to compute the distance from the current solution to g and ensure that the search moves towards the target.

In particular, consider routing problems, where a set of clients must be visited by a number of agents located at the depot along routes of minimum cost. We look for suggestions at the class of routing problems to which the m-TSP belongs, since, as mentioned, our application shares the constraints set with m-TSP. The most common move on which most neighborhoods are based is the relocate-move, which moves a client from a route to another. In fact, this move is able to transform any solution into another in a limited number of iterations. However, to use this move in a PR framework one should establish a correspondence between the routes in g and those in r based on an arbitrary naming convention, in order to identify which attributes of r are either identical to or different from g and relocate those clients currently served by a route to which they do not belong in g.

This issue has impaired a widespread application of PR techniques to routing problems: in particular, the only applications of PR to the classical Vehicle Routing Problem we are aware of are [23] and [42]; other studies related to routing problems involving PR concern arc routing problems [36, 44], clustered TSP [29], team orienteering [43], machine scheduling problems [45], project scheduling [35] capacitated location-routing problems [33], and periodic multi-depot [34] just to name a few. More in detail, in [23] PR is integrated in a Tabu Search, showing advantages with respect to a pure Tabu Search. The relocate move is used and a route-to-route correspondence is established by solving a matching problem between the routes in g and those in r based on the number of common clients; also [42] adopts the relocate move; a relocate-based distance between two solutions is defined as the minimum number of relocation moves required to transform one solution into the other, and a restricted neighborhood is defined accordingly. The resulting PR procedure is integrated within a Greedy Randomized Adaptive Search Procedure hybridized with a Variable Neighborhood Search. However, an extensive computational campaign did not prove any improvement due to PR; more surprisingly, the authors report that the overall procedure would perform better if PR would return the solution in the middle of the path between r and g rather than the best. This study is quite discouraging regarding the use of relocate moves or related ones that rely upon route-to-route correspondence. Indeed, our neighborhood based on the solution topology gets rid of this requirement and is based on a smaller piece of information, that is which node is the predecessor of a given node, as detailed in 3.4. This is closely tied to the move used in [44], however we developed this strategy independently based on i) an alternative topology based representation of the solution in the GA studied in [12], ii) a tentative to capture a feature of the application according to which the effects of the activation of a device are mostly influenced by the devices operated right before it.

3.3 Selection of reference candidates

As stated before, the reference set can be built up from the final populations of the GAs. In particular, diversity from the target beside quality should be taken into account, since the number of inspected solutions grows with the distance to the target; thus, different metrics can be combined together to properly filter the initial dataset and provide a set of good quality solutions not too close to each other.

The distance between two solutions can be evaluated considering either the routes of the teams or the activation times. In former studies diversity is measured on the graph representation of the routing problems, as discussed in 3.2. On the contrary, in this application we opt for measuring the diversity over the time representation: the graph representations of the solution would introduce a huge amount of redundancy [7], so the same vector of activation times can be mapped into different trees that may differ a lot with respect to the metrics defined for graph representations. On the opposite, in this problem, what makes a solution different from another one is the difference between the associated schedule, since the objective function depends on this. There are a few choices to measure the difference between two vectors; in particular, the metrics here proposed for the time representation are the Hamming distance $h (g, r) = \sum_{i = 1}^{N_{dev}} d_{i}$ (2) where g_i (r_i) is the activation time of device i in solution g (r) and d_i = 1 if g_i ≠ r_i while d_i = 0 otherwise, and the euclidean distance $e (g, r) = \sqrt{\sum_{i = 1}^{N_{dev}} (g_{i} - r_{i})^{2}}$ (3) between two vectors of activation times g and r. The former gives a measure about how many elements differ in the vectors, whereas the latter measures how much the vectors differ in terms of activation times. In order to prevent the inclusion of too similar vectors in rs, two thresholds β and γ are defined: given the target g a solution r is included only if h (g, r) ≥ β and e (g, r) ≥ γ.

As far as the quality of the reference set is concerned, a proper metric is to consider only solutions having quality within a certain percentage distance δ from the target’s one, i.e., $\frac{q (r) - q (g)}{q (g)} \leq δ$ holds for any r ∈ rs.

Finally, the choice of β, γ, and δ is really important, even more whenever the number of evaluations is limited. In fact, in this case excluding a promising solution may significantly affect the effectiveness of the approach.

3.4 A topology based path relinking

Building on the structure of the feasible region, i.e., the one of the mTSP, the few contributions in the literature about path relinking for routing problems may provide useful suggestions [42]. Path relinking for routing problems works on topological representations of routes, in which any route is expressed by an ordered set of visited customers [23 , 34]. For 3 teams v₁,v₂, and v₃, and 7 customers, namely c₁, …, c₇, a feasible solution assigns a route to each vehicle, e.g., v₁ = {c₁, c₂}, v₂ = {c₃, c₄}, v₃ = {c₅, c₆, c₇}. Equal solutions visit the customers along the same routes. To transform a solution into a different one, every customer should be relocated into the right position of the right route. In path relinking for routing problems, this is done iteratively by relocating one customer at each step.

The same representation can be used for mTSP. Since the routes have no names, the order of the devices within the routes is the only valuable information; moreover, in a route any device follows either the depot or a unique device, thus the order of activation within the routes can be stated by listing the devices’ predecessors, i.e., a list of tuples “(h_p, h_s)” meaning that device h_p precedes device h_s in the solution. For example, in the solutions in Fig. 2, three teams work overall on seven hydraulic devices, namely 1, …, 7; the initial solution r and the guiding solution g can be represented by the following predecessor lists:

$\begin{matrix} P_{r} = {(d, 3), (d, 4), (d, 6), (1, 5), (4, 7), (6, 1), (7, 2)}, \\ P_{g} = {(d, 3), (d, 4), (d, 6), (1, 7), (3, 5), (4, 1), (5, 2)} . \end{matrix}$ In general, two solutions a and b are equal iff P_a = P_b; whereas, whenever two solutions differ, the predecessors of a that are not in b are P_a–b = P_a \ P_a ∩ P_b, vice versa for b is P_b–a = P_b \ P_a ∩ P_b. Moreover, the cardinality card (P_a–b) measures the distance of a from b, with card (P_a–b) = card (P_b–a). For example, for r and g we have that $P_{r - g} = {(1, 5), (4, 7), (6, 1), (7, 2)},$ and $P_{g - r} = {(1, 7), (3, 5), (4, 1), (5, 2)};$ so, the distance between r and g is 4.

To get closer to b starting from the configuration of a, at least one device h_s such that (h_p, h_s) ∈ P_b–a should be relocated after its predecessor h_p in b; in this sense, the set P_b–a contains the possible moves to transform a into b. For example, (5, 2) ∈ P_g–r means that device 2 needs to be relocated right after 5, making r more similar and one step closer to g. This means that the neighbourhood of r with respect to g is the set of solutions $N (r, g) = {r_{n} ∣ c a r d (P_{g - r} \ P_{g - r_{n}}) = 1}$ . In other words, the neighbourhood of r with respect to g is the set of solutions r_n obtained by relocating one device in r according with P_g–r.

At the kth iteration, PR has to choose which device $h_{s} : (h_{s}, h_{p}) \in P_{g - r}^{k}$ is relocated, in order to move to r^k+1. To make this choice it evaluates by EPANET every $r_{n}^{k} \in N^{k}$ , and moves to the one with the lowest volume; this solution becomes the reference of the next step, and it is called r^k+1.

Every time a device has been relocated after its new predecessor their link becomes permanent and no further move can break it. Thus, whenever a device is chosen to be relocated after another one, it carries the following chain of fixed edges along with it; this prevents the current choice to destroy the previous ones. To implement this behaviour the procedure should be enriched with a memory, which stores the previous moves.

Figure 2 shows a possible path from r to g consisting of 4 intermediate steps. Table 1 reports at any step the values of P_g-r, the chosen move to r^k+1, and the fixed edges, for the example in Fig.2. The first move transforms the initial solution into r¹ by relocating 2 after 5; this edge is now fixed (shown in bold in the figure) and this move is stored into the memory, represented by a the dashed box. The second move from r¹ to r² relocates the chain 5 – 2 after 3. Then 7 is relocated after 1 achieving r³. Finally 1 is relocated after 4 together with its fixed successor 7. The target is finally reached.

Recall that in this real application only solutions with 3 routes are considered to be feasible since only 3 teams are available. So we exclude from P^k the moves that vary the number of routes, i.e., moves that either empty a route or add a new route.

This version of PR moves at most N_dev times and calls EPANET at most $\frac{N_{dev} (N_{dev} + 1)}{2}$ times. Sometimes the procedure visits the same solution twice or more, in such a case the solution would be evaluated by EPANET more times, wasting precious computing resources. For this reason the solving architecture is enriched with a cache, which stores the explored solutions and allows for saving a call to EPANET. It is worth noting that the procedure may explore the entire path between initial and guiding before the maximum number of EPANET calls expires. In such a case, the procedure selects another reference, and iterates over it. The algorithm continues until it reaches the maximum number of EPANET calls or reference solutions.

From now on, we refer to this version as the “routing” PR (PRseqr).

3.5 The time based variant

Another representation for the mTSP encodes a solution as a vector of activation times [13]. Given the feasible solutions r and g, the indexes of the differing elements is given by I_r–g = {i ∣ r_i ≠ g_i}. If r equals g then I_r–g =∅. To transform r into g iteratively, at each step k one element in I_r–g should be fixed to its value in g. This decreases by one the distance between r^k and g. Let r^k be the reference vector at the k-th step, I^k = I_r^k–g, and let F^k = I_r–g \ I^k be the set of indexes that have been already fixed. The next solution r^k+1 in the path between r and g is obtained by keeping $r_{f}^{k + 1} = g_{f}$ for all f ∈ F^k and fixing the new element $r_{i}^{k + 1} = g_{i}$ for one i ∈ I^k.

If the remaining elements of r^k+1 were the same as in r (or r^k) the resulting vector could not correspond to a feasible schedule. So, these elements are allowed to change and are chosen by solving a constrained optimization problem whose constraints define a mTSP, where the elements in F^k ∪ {i} are fixed whereas the other (non-fixed) elements are the actual integer variables of the model. The objective is to optimize these variables, so that their values are as close as possible to the ones in r. To do that, the model minimizes the Euclidean distance of the non-fixed elements from r in a manner similar to what it has been been done in [13], i.e., given i ∈ I^k: $dist (r^{k}, r) = \sum_{j \in I^{k} \ {i}} | r_{j}^{k + 1} - r_{j} |$ (4)

Notice that a feasible vector always exists, being g a feasible solution of the program. The neighbourhood of r^k is then defined as follows: $N^{k} = {r^{k + 1} ∣ c a r d (I^{k} \ I^{k + 1}) = 1, minimizes (4)} .$

The procedure explores every solution by varying the index i ∈ I^k, and for each i it calls the optimiser to compute a new feasible vector, finally it evaluates the solution by calling EPANET. The solution in N^k having the lowest volume is selected to be the reference solution for the next step. Figure 3 shows, on a graph with 7 devices, how routes change when an additional element becomes fixed; e.g., in r⁶ the activation time “(26)”, which was (20) in r⁵, has been fixed, and the related device is now visited by another route. The minimization of dist (r⁶, r) also transforms (28) in r⁵ into (29) in r⁶, by changing the route of the concerned device; this new activation time is clearly very close to (27) in r.

The constrained optimization program minimizing(4) can be stated by any declarative paradigm, such as Constraint Programming [10], Constraint Logic Programming [24], Answer Set Programming [16, 26], Mixed Integer Linear Programming [28]; so some suitable solvers are: Gecode [15], ECLiPSe [40], DLV [27], Clasp [14], SCIP [1], Gurobi [21].

Since the number of feasible moves decreases at each iteration of at least one unit, the maximum number of EPANET calls is again $\frac{N_{dev} (N_{dev} + 1)}{2}$ . Also this PR uses a cache to store the explored solution, so it ends up whenever either no more EPANET calls are available, or rs is empty.

Notice how this technique integrates MILP, hydr-aulic simulators, and PR into the solving architecture; thus, we will refer to it as the “hybrid” PR (PRh).

4 Computational results

We extended the experimental campaign introduced in [31] keeping the same experimental platform. In the previous work hybrid path relinking was in general better performing than the route based one in the 20 contamination scenario tested; however, a deeper insight about this result can be obtained by extending the experiments on more scenarios.

The experiments were performed on the Ferrara’s hydraulic network, which supplies drinking water to about 130, 000 inhabitants. 42 contamination scenarios (1 … 42) were tested, 22 more than the previous work [31]. Basing on the techniques proposed in [20], 3 teams of technicians were considered to be available to operate on 13 hydraulic devices, namely 7 valves and 6 hydrants. The hydraulic simulator we used was EPANET [39], and it takes about 5 seconds to evaluate a schedule of the selected devices on each contamination scenario. Even though EPANET is open-source, the simulation procedures and the network specifications are sensitive data for hydraulic engineers and can not be disclosed. Due to the running time required by a simulator call, the results here presented requested about 32 days of computation on AMD Athlon 64 X2 Dual Core Processor 5600+ and 2800MHz.

Genetic algorithms proposed in [13] were allowed a maximum of 500 EPANET calls, and the population was sized to 20 individuals. These values were calibrated in previous works, and the GAs typically converge within the 500 EPANET calls. As mentioned, we observed that the final solution depends on the initial population as typically GA gets stuck at a local optimum, so parallel small sized independent GAs explore the search space better than one big sized GA and can be parallelized. In this study, Path Relinking is tested to explore the region enclosing such solutions.

We experimentally verified that with 5 GAs the reference set in output was less diverse than with 10 GAs; in fact, PR performed better in the second case. This complies with the fact that PR is more effective when the reference set has high quality and diversity (see Section 3).

Regarding the filtering of the reference set out of the union of the final populations of the 10 GAs, preliminary experiments led us to the following values, which provided good results for all the scenarios: β = 11, γ = 1, and δ = 2. This is a very light filtering that guarantees a sufficiently large reference set; the success of this policy may be explained by the following facts: the selection of the initial and the guiding solutions is driven by the objective function value (the lower the better), so that the limited number of solution evaluations allowed leaves out most of the worst solutions; the GA population size coincides with the size of the reference set suggested by Glover [18] (20), and the final populations of the 10 GAs have several elements in common; even short paths are worth being explored since the objective function is sensitive to small variations of the schedule; finally, even though the reference set may contain solutions very similar to each other and this may lead to search trajectories that share some points this is not penalizing due to the cache mechanism. It is worth noting that in our procedure the filtering is applied once for all, at the beginning, and the reference set is never updated during the procedure. On the contrary, when the reference set is updated with the new solutions, dynamic and adaptative filtering can be implemented in order to preserve the diversity and the cardinality of the reference set.

The hypotheses tested in the experimental campaign are whether:

intensificating over the final population of 10 GAs may improve the local optima;

an intensification step may be more effective than an additional independent GA run;

given that, hybrid time based GAs outperformed route based ones in [13], likewise we expect (and preliminary proved in [31]) that time based path relinking will dominate the route based one.

To tests such hypotheses the path relinking procedures were started from the final populations of 10 independent GAs. The two PRs were run independently compared to an additional independent GA run, this to be considered a strengthening run of the same first 10. In this way all the approaches are directly comparable in terms of computational resources. The resulting experimental platform is depicted in Fig. 4. It shows that the 10 GAs’ final populations are merged and filtered to yield the reference set rs; PRr and PRh are started from that rs and return the contaminated water volumes V^r and V^h, respectively; V^g is the best over the solution returned by the additional GA and the former 10 (so the best in rs): basically it is like running 11 independentGAs.

To weaken the randomness, the tests were repeated 10 times on each contamination scenario. To disambiguate, the 10 independent GA runs are intended to be local, with respect to the 10 general runs that are called global. PRr, PRh, and the additional GA were allowed 500 EPANET calls each; the total number of simulations is then 6500 for any configuration. Consequently, the total number of EPANET calls is 6500 · 42 · 10 = 2, 730, 000, corresponding to about 32 days of (a single) CPU time. The constrained optimization model minimizing (4) was implemented as a Mixed Integer Linear Programming and solved by Gurobi [21]; its solving time is negligible within the single run (see [13] for some details). Increasing the number of devices to be scheduled would increase the computing time either; however, this number can be minimized [20] and the resulting instance is solved by the MILP solver within a reasonabletime.

We introduce the following benchmark metrics for the first optimization stage:

$s_{c}^{*} = m i n {r s_{c}}$ : the best of the c-th trial of 10 GAs;

$s * = m i n_{c} {s_{c}^{*}}$ : global optimum within the 10 × 10 GA runs.

Comparing the path relinking outcome with the former assesses whether and how much the intensification step has been effective within the single trial; comparing with the latter instead assesses how effective is the intensification step in finding very good solution with respect to the GAs’ ones. To make this comparison we introduce the following metrics for each method of the second optimization stage:

V_c: the contaminated water in output from one of the three method in the c-th trial;

$Δ_{c} = s_{c}^{*} - V_{c}$ : how much contaminated water is saved by the second step with respect to the reference set: as $V_{c} \leq s_{c}^{*}$ , Δ_c ≥ 0 always holds;

$\bar{Δ}$ : the averaged amount over the trials.

These quantities are measured in litres (l). Other useful metrics are:

Ic= {1 ⇔ Δ_c > 0, 0 otherwise}: whether the method improves the local optimum;

ΣI = ∑_cI_c: how frequently the method can improve the starting solution;

I^*= ∑_c {1 ⇔ min {V_c} < s^*, 0 otherwise}: how many times the method improves the global optimum (overall the 10 trials).

To better explain these metrics Table 2 reports the values for the scenario 31 for each trial c. In this scenario the additional GA never improved s^*, PRr improved it four times yielding the best $\bar{Δ}$ , whereas PRh improved it seven times and found the best solution in the whole data result. So it is not infrequent that a method may compute the best solution in the 10 trials while another method may yield better $\bar{Δ}$ . In fact PRr was able to improve $s_{c}^{*}$ with high values even though it did it less times than PRh. PRh improved s^* once in run 9.

For the sake of readability, we aggregate the results by scenarios in Table 3 as done in rows h-r in Table 2. To test the first hypothesis, recall we compare how many times the three approaches improve $s_{c}^{*}$ , and how many litres they save, that are outcomes of the first optimization step. The additional GA (+1 GA) improved $s_{c}^{*}$ only 9 times, but on 9 different scenarios. PRr improved $s_{c}^{*}$ overall 54 times and in 21 scenarios never improved it. PRh improved $s_{c}^{*}$ overall 272 times, at least twice in every scenario. So, on average PRr improved $s_{c}^{*}$ 1.2 times of 10, while PRh did it 6.5 times; also, that one time PRr improved the local optimum it saved, on average, 20.9 litres of contaminated water, while PRh saved 63.8 litres. This suggests that intensificating pays back, and does it more than an additional (non-intensificating) run of GA, which proves also the second hypothesis.

Notice that one could consider to run either the new 10 (additional) GA or the PRr or the PRh in the same parallel architecture as the first 10 GAs were, considering them as an entire second stage optimization process. In such a case one would select at the end the best schedule min {V_c} in order to compare it with s^*. By looking at the component (I^*) it turns out that the second stage of GAs improved s^* just once in scenario 17 (which were pretty unlucky); PRr improved s^* in 5 scenarios, and PRh improved it in 26 scenarios; in 14 scenarios out of 42 neither PRr nor PRh improved. Also, the last row shows that PRh is expected to improve s^* about once in any scenario. Consequently, in this context intensification by path relinking is better performing than rerun the first stage again.

As far as the third hypothesis to be tested, we compare in how many scenarios the PRr improved s^* (I^*) whereas PRh did not, which happened in scenario 1; on the contrary, in 23 scenarios PRh improved s^* whereas PRr did not. We can also compare for how many scenarios PRr got better values of I than PRh and vice versa, obtaining that only in scenario 1 PRr improved the local optimum more times than PRh (8 vs 7). Moreover, PRh on average improved $s_{c}^{*}$ six as many times as PRr (see last row). In terms of contaminated water, PRh saved 3 times the volume than PRr (63.8 l vs 20.9 l). This suggests that PRh more likely finds better schedules than PRr, and decreases the contaminated water more than PRr. Despite of the fact that PRh is generally better performing than PRr, there is not a total dominance, which suggests further investigations.

Figures 6 –8 show three boxplots for each contamination scenario. The legend of these charts is provided in Fig. 5 the x-axis reports the amount of saved contaminated water Δ_c, while along the y-axis three horizontal boxplots are listed for any scenario: the GA one, then PRr, and finally PRh, from the top to the bottom. Outliers may appear as points.

These charts show many outliers in many scenarios. This means that path relinking has some chances to find particularly good schedules; in other words, these schedules are within the spectrum of the solutions explored by the path relinking but cannot be reached. This suggests to empower the path relinking to use the allowed 500 evaluation more wisely.

5 Conclusions and future works

This work extends the experimental analysis of an intensification step to be performed by path relinking on the local optima in output by 10 parallel genetic algorithms. We tested 42 contamination scenarios, 22 more than the previous work [31] every test was repeated 10 times to deal with randomness, and the total computing time was about 32 days in a single CPU due to the great deal of effort in hydraulic simulations that was needed to compute the schedules’ quality. Results show that intensificating pays back. In fact path relinking is able to decrease the amount of contaminated water and it achieves this result more often than an additional independent (thus non-intensificating) run of a genetic algorithm. The time based path relinking, which integrates a m-TSP solver, outperforms route based path relinking in many scenarios: it improves the local optimum more often and saves more contaminated water; however, there is no clear dominance of time based over the route based and further investigations should not neglet the latter. Indeed, in some specific scenario the route based path relinking finds very good solutions not matched by those reached by the time based one. The outliers in the boxplots suggests that the search capability of path relinking can be further improved to perform the exploration more in breadth, given the same number of solution evaluation: more extensively, a smarter strategy should be devised to exploit the available knowledge accumulated in the current as well as in the previous paths to prioritize some search trajectories other than others, so to reach the highest quality solutions known to exists within the allotted 500 trials. Moreover, well-known improved variants of path relinking can be applied [37]: one idea is to interrupt the search when the previous moves were not good, implementing the so called truncated path relinking, or the path direction could be reversed swapping the initial and guiding solutions. These hints plus a comparison with intensificating genetic algorithms are currently under study.

Footnotes

Acknowledgments

We thank Stefano Alvisi and Marco Franchini for assistance with the hydraulic simulator and the instance they developed, and for the fruitful discussions about the hydraulic engineering.

References

1.
Achterberg T. , SCIP: Solving constraint integer programs, Mathematical Programming Computation, 1(1) (2009), 1–41. http://mpc.zib.de/index.php/MPC/article/view/4.

2.
Alvisi S. , Franchini M. , Gavanelli M. , Nonato M. , Nearoptimal scheduling of device activation in water distribution systems to reduce the impact of a contamination event, Journal of Hydroinformatics, 14(2) (2012), 345–365.

3.
April J. , Glover F. , Kelly J.P. , Laguna M. , Simulationbased optimization: Practical introduction to simulation optimization, In Proceedings of the 35th Conference on Winter Simulation: Driving Innovation, WSC ’03, pages 71–78. Winter Simulation Conference, 2003.

4.
Bektas T. , The multiple traveling salesman problem: an overview of formulations and solution procedures, Omega, 34(3) (2006), 209–219.

5.
Bent R. , Coffrin C. , Judi D. , McPherson T. and van Hentenryck P. , Water distribution expansion planning with decomposition, In WDSA 2012:14th Water Distribution Systems Analy-sis Conference, 24-27 September 2012 in Adelaide, South Aus-tralia, page 305. Engineers Australia, 2012.

6.
Blum C. and Roli A. , Metaheuristics in combinatorial optimization: Overview and conceptual comparison, ACM Com-puting Surveys (CSUR), 35(3) (2003), 268–308.

7.
Carter A.E. and Ragsdale C.T. , A new approach to solving the multiple traveling salesperson problem using genetic algorithms, European Journal of Operational Research, 175(1) (2006), 246–257.

8.
Cattafi M. , Gavanelli M. , Nonato M. , Alvisi S. and Franchini M. , Optimal placement of valves in a water distribution network with CLP(FD), Theory and Practice of Logic Programming, 11(4-5) (2011), 731–747.

9.
Feo T.A. and Resende M. G. , Greedy randomized adaptive search procedures, Journal of global optimization, 6(2) (1995), 109–133.

10.
Frühwirth T. and Abdennadher S. , Essentials of Constraint Programming. Springer, 2003.

11.
Gavanelli M. , Nonato M. and Peano A. , An ASP approach for the valves positioning optimization in a water distribution system, Journal of Logic and Computation, 25(6) (2015), 1351–1369.

12.
Gavanelli M. , Nonato M. , Peano A. , Alvisi S. and Franchini M. , Genetic algorithms for scheduling devices operation in a water distribution system in response to contamination events, In Hao J.-K. , and Middendorf M. , editors, Evolutionary Computation in Combinatorial Optimization, volume 7245 of Lecture Notes in Computer Science, pages 124–135. Springer Berlin/Heidelberg, 2012.

13.
Gavanelli M. , Nonato M. , Peano A. , Alvisi S. and Franchini M. , Scheduling countermeasures to contamination events by genetic algorithms, AI Communications, 28(2) (2015), 259–282.

14.
Gebser M. , Kaufmann B. , Kaminski R. , Ostrowski M. , Schaub T. and Schneider M. , Potassco: The Potsdam Answer Set Solving Collection, AI Communications, 24(2) (2011), 107–124.

15.
Gecode Team, Gecode: Generic constraint development environment, 2006 Available from http://www.gecode.org.

16.
Gelfond M. , Answer sets, In van Harmelen F. , Lifschitz V. , and Porter B. , editors, Handbook of Knowledge Representa-tion chapter 7, pages 285–316. Elsevier Science, 2008.

17.
Glover F. , Heuristics for integer programming using surrogate constraints, Decision Sciences, 8(1) (1977), 156–166.

18.
Glover F. , Laguna M. and Martí R. , Fundamentals of scatter search and path relinking, Control and Cybernetics, 39 (2000), 653–684.

19.
Gonçalves J.F. , de Magalhães Mendes J.J. and Resende M.G.C. , A hybrid genetic algorithm for the job shop scheduling problem, European Journal of Operational Research, 167(1) (2005), 77–95.

20.
Guidorzi M. , Franchini M. and Alvisi S. , A multi-objective approach for detecting and responding to accidental and intentional contamination events in water distribution systems, Urban Water, 6(2) (2009), 115–135.

21.
Gurobi Optimization, Inc, Gurobi optimizer reference manual, 2014. http://www.gurobi.com.

22.
Haq Z.U. and Anwar A.A. , Irrigation scheduling with genetic algorithms, Journal of Irrigation and Drainage Engineering, 136(10) (2010), 704–714.

23.
Ho S. and Gendreau M. , Path relinking for the vehicle routing problem, Journal of Heuristics, 12(1-2) (2006), 55–72.

24.
Jaffar J. and Maher M.J. , Constraint logic programming: A survey, Journal of Logic Programming, 19/20 (1994), 503–581.

25.
Krause A. , Leskovec J. , Guestrin C. , VanBriesen J. and Faloutsos C. , Efficient sensor placement optimization for securing large water distribution networks, Journal of Water Re-sources Planning and Management, 134(6) (2008), 516–526.

26.
Leone N. , Logic programming and nonmonotonic reasoning: From theory to systems and applications, In Baral C. , Brewka G. , and Schlipf J. , editors, Logic Programming and Nonmonotonic Reasoning, volume 4483 of Lecture Notes in Computer Science, pages 1–1. Springer Berlin Heidelberg, 2007.

27.
Leone N. , Pfeifer G. , Faber W. , Eiter T. , Gottlob G. , Perri S. and Scarcello F. , The DLV system for knowledge representation and reasoning, ACM Transactions on Computational Logic (TOCL), 7(3) (2006), 499–562.

28.
Martin R. , Large Scale Linear and Integer Optimization: A Unified Approach. Springer US, 1999.

29.
Mestria M. , Ochi L.S. and de Lima Martins S. , GRASP with path relinking for the symmetric euclidean clustered traveling salesman problem, Operations Research, 40(12) (2013), 3218–3229.

30.
Murray R. , Hart W. , Phillips C. , Berry J. , Boman E. , Carr R. , Riesen L.A. , Watson J.-P. , Haxton T. , Herrmann J. , Janke R. , Gray G. , Taxon T. , Uber J. and Morley K. , US environmental protection agency uses operations research to reduce contamination risks in drinking water, Interfaces, 39(1) (2009), 57–68.

31.
Nonato M. and Peano A. , Path relinking for a constrained simulation-optimization team scheduling problem arising in hydroinformatics, In Gavanelli M. , Lamma E. and Riguzzi F. , editors, AI*IA 2015 Advances in Artificial Intelligence, volume 9336 of Lecture Notes in Computer Science, pages 31–44. Springer International Publishing, 2015.

32.
Peano A. , Solving Real-Life Hydroinformatics Problems with Operations Research and Artificial Intelligence. PhD thesis, University of Ferrara, 2015.

33.
Prins C. , Prodhon C. and Wolfler R. , Calvo, Solving the capacitated location-routing problem by a grasp complemented by a learning process and a path relinking, 4OR, 4(3) (2006), 221–238.

34.
Rahimi-Vahed A. , Crainic T. , Gendreau M. and Rei W. , A path relinking algorithm for a multi-depot periodic vehicle routing problem, Journal of Heuristics, 19(3) (2013), 497–524.

35.
Ranjbar M. , Kianfar F. and Shadrokh S. , Solving the resource availability cost problem in project scheduling by path relinking and genetic algorithm, Applied Mathematics and Computation, 196(2) (2008), 879–888.

36.
Reghioui M. , Prins C. and Labadi N. , Grasp with path relinking for the capacitated arc routing problem with time windows, In Giacobini M. , editor, Applications of Evolutionary Comput-ing, volume 4448 of Lecture Notes in Computer Science pages 722–731. Springer Berlin Heidelberg, 2007.

37.
Ribeiro C.C. and Resende M.G. , Path-relinking intensification methods for stochastic local search algorithms, Journal of Heuristics, 18(2) (2012), 193–214.

38.
Ribeiro C.C. and Vianna D.S. , A hybrid genetic algorithm for the phylogeny problem using path-relinking as a progressive crossover strategy, International Transactions in Operational Research, 16(5) (2009), 641–657.

39.
Rossman L.A. , EPANET 2 users manual, National Risk Management Research Laboratory, Office of research and development, U.S. Environmental Protection Agency, USA., 2000.

40.
Schimpf J. and Shen K. , Eclipse - from LP to CLP, TPLP, 12(1-2) (2012), 127–156.

41.
Simonis H. , Constraint applications in networks, Handbook of constraint programming, 2 (2006), 875–903.

42.
Sörensen K. and Schittekat P. , Statistical analysis of distance based path relinking for the capacitated vehicle routing problem, &}, Operations Research, 40(12) (2013), 3197–3205.

43.
Souffriau W. , Vansteenwegen P. , Berghe G.V. and Oudheusden D.V. , A path relinking approach for the team orienteering problem, Computers & Operations Research, 37(11) (2010), 1853–1859. Metaheuristics for Logistics and Vehicle Routing.

44.
Usberti F.L. , França P.M. and França A.L.M. , GRASP with evolutionary path-relinking for the capacitated arc routing problem, Computers & Operations Research, 40(12) (2013), 3206–3217.

45.
Vallada E. and Ruiz R. , Genetic algorithms with path relinking for the minimum tardiness permutation flowshop problem, Omega, 38(1-2) (2010), 57–67.