A matheuristic approach for the Pollution-Routing Problem

This paper deals with the Pollution-Routing Problem (PRP), a Vehicle Routing Problem (VRP) with environmental considerations, recently introduced in the literature by [Bektas and Laporte (2011), Transport. Res. B-Meth. 45 (8), 1232-1250]. The objective is to minimize operational and environmental costs while respecting capacity constraints and service time windows. Costs are based on driver wages and fuel consumption, which depends on many factors, such as travel distance and vehicle load. The vehicle speeds are considered as decision variables. They complement routing decisions, impacting the total cost, the travel time between locations, and thus the set of feasible routes. We propose a method which combines a local search-based metaheuristic with an integer programming approach over a set covering formulation and a recursive speed-optimization algorithm. This hybridization enables to integrate more tightly route and speed decisions. Moreover, two other"green"VRP variants, the Fuel Consumption VRP (FCVRP) and the Energy Minimizing VRP (EMVRP), are addressed. The proposed method compares very favorably with previous algorithms from the literature and many new improved solutions are reported.

Introduction
The Vehicle Routing Problem (VRP) and its variants have been the subject of considerable research efforts in the past year, mostly due to the growing number of additional constraints and objectives arising from real-world problems (Vidal et al., 2013a). Given the growing global concern about environmental issues, VRPs have recently started to incorporate "green" aspects such as pollution and alternative fuels, among others (see, e.g., the recent review of Lin et al. 2014 on VRPs with environmental issues).
Environmental costs due to greenhouse gas (GHG) emissions are not usually paid directly by the companies. Nevertheless, some countries are developing emissions trading schemes (e.g., the European Union Emissions Trading System) to make the companies responsible for product between load and distance for each arc. Similar approaches, i.e, those that make use of the vehicle load to minimize the fuel consumption or CO 2 emissions, were presented by Peng e Wang (2009), Scott et al. (2010), Ubeda et al. (2011) and Xiao et al. (2012). The latter introduced the Fuel Consumption Vehicle Routing Problem (FCVRP). Kara et al. (2008) presented a VRP with cumulative costs, which generalizes the EMVRP as well as the Minimum Latency Problem and the m-Traveling Repairman Problem. Kopfer et al. (2013) took the load into account for CO 2 emissions evaluations in the presence of heterogeneous vehicle types.
Minimizing the fuel consumption considering only the load and distance can be insufficient since the travel speed plays a major role. This speed is directly affected by the road congestion. In view of this, Kuo (2010) proposed a model for minimizing the total fuel consumption where the speeds are time-dependent and the load is used to estimate the cost. A Simulated Annealing algorithm was implemented to solve a set of instances from the Solomon benchmark (Solomon, 1987). Later, Kuo e Wang (2011) devised a Tabu Search algorithm for the same problem. Other time-dependent problems with emission minimization can be found in Figliozzi (2011), Saberi e Verbas (2012) and Jabali et al. (2012), and many other references relevant to green logistics are mentioned in the surveys of Dekker et al. (2012), Salimifard et al. (2012), Lin et al. (2014) and Demir et al. (2014b). Bektaş e Laporte (2011) proposed the Pollution-Routing Problem (PRP), which seeks to minimize both operational and environmental costs, taking into account the customers timewindow constraints. The total travel distance, the amount of load carried per distance unit, the vehicle speeds and the duration of the routes are the main cost components. Three different variants were also presented, considering either distance, weighted load and energy minimization. The authors performed an extensive experimental analysisto capture the trade-off between each variation, as well as the effect of the speed and time-window constraints on the distance, energy and costs.
The PRP was addressed with a two-phase heuristic in Demir et al. (2012). In the first phase, the VRPTW is solved by means of an Adaptive Large Neighborhood Search (ALNS), including five insertion operators and twelve removal operators. In a second phase, vehicle speeds are optimized using a recursive algorithm. Computational experiments were carried out for instances with up to 200 customers. Other recent developments consider generalizations of the problem. A bi-objective variant considering fuel and driving time minimization is presented in Demir et al. (2014a), andFranceschetti et al. (2013) consider the time-dependent PRP.
It should be noted that several aspects of the PRP are conflicting. Higher speeds imply routes with shorter durations, but at the same time result in a larger amount of emissions and vice-versa. Hence, to reduce pollution, speed on arcs may be decreased to be closer to the speed which minimizes emissions. Yet with lower speed, the set of feasible VRP routes may be empty or drastically smaller. As a consequence, an optimal solution of the reduced speed VRP can have longer distance and even more emissions in some cases. One main goal of our study was to better integrate route and speed decisions in order to find quickly a suitable balance between these antagonist aspects.

Problem description
The PRP can be defined as follows. Let G = (V, A) be a complete and directed graph with a set V = {0, 1, 2, . . . , n} of vertices and a set A = {(i, j) ∈ V 2 , i = j} of arcs. Vertex 0 represents the depot where a fleet of m identical vehicles with capacity Q is based. Vertices V − {0} correspond to customers, characterized by a non-negative demand q i for a single product, a service time τ i and a specified time-window interval [a i , b i ] for service. We assume that q 0 = 0 and τ 0 = 0. Each arc (i, j) ∈ A represents a travel possibility from node i to j for a distance d ij .
A particularity of the PRP is that the speed v ij on each (i, j) is itself a decision variable, valued between [v min , v max ]. Indeed, each vehicle emits on a certain amount of GHG which depends of weight and speed, among other factors. The PRP aims to find a speed matrix (v) ij and a set of routes R (such that |R| ≤ m) to serve all customers while minimizing environmental and operational costs. Each route σ ∈ R, σ = (σ 1 , . . . , σ |σ| ) starts and ends at the depot, i.e., σ 1 = 0 and σ |σ| = 0, the total demand on each route should not exceed the vehicle capacity, and every customer must be visited within its time window.
Equation (1) defines for a route σ the vehicle load f σ i σ i+1 when traveling on arc (σ i , σ i+1 ), and Equation (2) defines recursively the arrival time t i to customers, knowing that each route starts at time zero, that it is allowed to arrive early at a customer location and wait the start of its time window a i , but that a late arrival is not permitted.
The PRP objective is based on the assumption that CO 2 emissions are approximately proportional to fuel consumption. Using the comprehensive emissions model of Barth et al. (2005), Scora e Barth (2006) and Barth e Boriboonsomsin (2008), Bektaş e Laporte (2011) Parameters w 1 , w 2 , w 3 , w 4 are based on fuel properties, vehicle and network characteristics. Their meaning and value are discussed in Appendix 1, based on Bektaş e Laporte (2011) and Demir et al. (2012). Finally, the overall objective of the PRP, given in Equation (4), considers both the fuel consumption with a cost of ω fc pounds (£) per liter, and the driving costs at a rate of ω fd per unit of time. The The minimum v * f on this function, i.e., the speed value that minimizes fuel costs, is given in Equation (5).
Similarly, for any arc (σ i , σ i+1 ), assuming that there is no waiting time in the route after σ i , the travel cost F fd σ i σ i+1 (v σ i σ i+1 ) including driver wages is given in Equation (6). The speed value that minimizes fuel and driver costs is then expressed in Equation (7). Both values, v * f and v * fd , are independent of the arc under consideration.
The two other problems considered in this work, EMVRP and FCVRP, do not take into account time-window constraints and speed decisions. In practice, these problems can be seen as CVRPs with green-oriented objective functions. The objective of the EMVRP is based on a simplified emission model, which takes solely into account the distance and the load×distance factor (Kara et al., 2007). It is expressed in Equation (8), where ω represents the weight of a vehicle without cargo.
Finally, the FCVRP is based on a linear function of fuel consumption per distance, which considers the vehicle's no-load/curb weight, the load carried, and the fixed cost of the vehicle. The objective is expressed in Equation (9), where h is the vehicle fixed cost, ρ * and ρ 0 are the fuel consumption rate of the vehicle without any load and fully loaded, respectively.
Using the same values of Xiao et al. (2012), with h = 0, ω fc = 1, ρ * = 2, and ρ 0 = 1, the objective can be reformulated as in Equation (10). It is noteworthy that, as in the case of the EMVRP, the FCVRP involves a linear combination of the distance and the load×distance factor.
4 The proposed ILS-SP-SOA matheuristic The proposed algorithm, called ILS-SP-SOA, combines an iterated local search with speed optimization procedures and integer programming optimization over a set partitioning formulation. As a blend of metaheuristic and exact procedures, this category of method is usually called matheuristic (Maniezzo et al., 2009). In ILS-SP-SOA, the initial solution construction and local search, as well as three perturbation procedures aim at minimizing the cost by considering routing decisions without changing the speeds. This sub-problem can be seen as a VRPTW with the objective of Equation (4). We demonstrate in Section 4.1 how local-search moves can be evaluated in O(1) time for this setting. The search is complemented by a recursive speed-optimization procedure, described in Section 4.3, which is applied on each local minimum of the "routing" local search. New speed decisions are included in a dynamic speed matrix, which is in turn used in the subsequent VRPTW sub-problems. Finally, the routes associated to local minimums are stored in a pool, and used by an integer optimization procedure over a set partitioning formulation to generate possible better solutions composed of a different recombination of routes. In this process, each route may be associated to a different speed matrix, and thus this exact procedure adds an additional level of integration between route and speed decisions.
The outline of ILS-SP-SOA is presented in Algorithm 1. The method performs n r restarts (Lines 3-23) of a hybrid procedure combining ILS, SOA and SP. At each restart, the speed matrix v is first initialized with the maximum speed, v ij = v max for all (i, j) ∈ A. A solution S is obtained by applying local search and SOA on an initial solution (Line 8). This initial solution is generated using the modified cheapest insertion heuristic of Penna et al. (2013), considering the PRP objective and the time-window relaxation of Section 4.1. The speed of arcs associated to S is then updated in the matrix v.
Algorithm 1 ILS-SP-SOA(n r , n ils , n sp , n pool , T mip , seed) 1: S best-all ← ∅; f (S best-all ) ← ∞; {S best-all is the overall best solution} 2: P perm ← ∅; i r = 0; { P perm is the pool of permanent routes in the SP} 3: while i r < n r do 4: i r ← i r + 1;

5:
S best ← ∅; f (S best ) ← ∞; {S best is the best solution of the restart phase} i ils ← i ils + 1;

21:
if f (S best ) < f (S best-all ) then

24:
if number of consecutive restarts ≥ n pool then 25: The ILS is then run until n ils successive iterations of local search and perturbation are reached without improvement. At each iteration (Lines 10-18) the local optimal solution is modified by one of the three perturbation mechanisms, selected at random with different probabilities as described in Section 4.2. Time-window constraints are relaxed when performing a perturbation move. This modified solution is possibly improved by applying local search and SOA (Line 12), and the speed matrix is updated (Line 13). A temporary pool of routes is updated during the ILS loop by adding routes associated to local optimal solutions (Line 14). If the number of ILS iterations without improvement is equal to n, the perturbation Change Speeds is applied, as explained later in Subsection 4.2. Finally, the speed matrix v is reinitialized when the number of ILS iterations without improvement is greater than n ils /2. (Lines 17-18).
If the size of the instance is smaller than a given parameter n sp , then the SP method is called only after the last restart phase; otherwise, it is called after every restart. The SP method attempts to create a new solution from routes of a temporary pool of routes P temp derived from local optimums of the local search, which is cleared after each n pool restarts (Lines 24-25), and from a permanent pool of routes P perm which contains the routes associated to the best solutions S best of each restart phase (Line 23). The SP problems are solved using a Mixed Integer Programming (MIP) solver which calls the ILS scheme every time a new incumbent solution is found (Line 20). This collaborative approach, described in Subramanian et al. (2013), not only potentially speeds up the solver runtime, but also may contribute to find better solutions than the best possible combination of routes in the pools. If this happens, the pool P temp is also updated. The MIP solver is run several consecutive times until no improvement is found over S best , and a time limit T mip is imposed to each MIP execution to avoid any unpredictable excessive CPU time. Finally, the algorithm returns the best solution found among all restarts (Line 26).

Local Search with efficient move evaluation
The local search procedure is based on the Randomized Variable Neighborhood Descent (RVND) of Subramanian et al. (2010). It relies on five inter-route neighborhoods, namely: Shift(1,0), Shift(2,0), Swap(1,1), Swap(2,2), 2-opt * and five intra-route neighborhoods, namely: Reinsertion, Or-opt2, Or-opt3, Exchange and 2-opt. The neighborhood structures based on Shift operators move one or more consecutive customers from a route to another one. Those based on Swap consists of interchanging one or more customers from one route with one or more customers from another. Reinsertion and Or-opt follow the idea of Shift, but involving customers of a single route. Exchange is the intra-route version of Swap(1,1). Finally, 2-opt is the classical Traveling Salesman Problem (TSP) intra-route operator, whereas 2-opt * is the interroute version of 2-opt. A detailed description of these neighborhoods, as well as the Auxiliary Data Structures (ADSs) used to enhance the performance the local search can be found in , Penna et al. (2013) andVidal et al. (2013b).
A local search for the VRPTW may have difficulties to generate feasible solutions, possibly compromising the convergence towards good solutions. To circumvent this issue, penalized infeasible solutions w.r.t. time windows are considered, and the objective includes "time-warp" penalties as in (Vidal et al., 2013b). To compute the cost of new routes generated by the local search in amortized constant time, other ADSs based on subsequences concatenation were implemented. For any subsequence σ the algorithm stores and maintains: • minimum duration T (σ), • minimum time-warp use T W (σ), • earliest E(σ) and latest visit L(σ) to the first vertex allowing a schedule with minimum duration and minimum time-warp use, For a subsequenceσ involving a single customer, i, these ADSs are computed as follows: If the first node in the sequence is a depot, e.g., ifσ = 0, then E(σ) = L(σ) = 0. This prevents a delayed departure. The following equations enable then to derive ADSs for larger subsequences obtained by concatenation ⊕: The penalized cost of a route σ can be derived from these structures as shown in Equation (23) where ω tw is the penalty for one unit of time warp.

Perturbation Mechanisms
The three mechanisms applied during the perturbation phase are described as follows.
• Shift to End -Shifts one customer from one route to the end of other route, randomly with uniform distribution.
• Change Speeds -This perturbation modifies the speed matrix by changing the speeds associated to arcs of one random route in the current best solution S best with one random speed from the set {v * f , v * fd , v max }, with uniform distribution.
During the perturbation phase, Shift to End and Merge routes are randomly selected with different probabilities, 90% and 10%, respectively. Any perturbation leading to an infeasible solution w.r.t. capacity constraints is undone and a new perturbation is attempted. Finally, the additional Change Speeds perturbation is applied after n consecutive ILS iterations without improvement, e.g., when i ils = n.

Speed Optimization Algorithm
Speeds have a direct impact on customers' arrival time, i.e., in meeting their time windows. When routes are fixed, the PRP leads to a speed optimization problem which consists of finding the optimal speeds for each arc while respecting customers' time windows. This problem is solved by a recursive algorithm (Algorithm 2). This method is an adaptation of the RSA procedure of Norstad et al. (2011) and Hvattum et al. (2013), also considering drivers wages, possible waiting times, and a non-fixed arrival time at the last customer.
Algorithm 2 is applied on the complete route, by setting s = 1 and e = |σ|. Let t σ i be the time to start the service of customer σ i . Early and late arrivals are not considered in the course of the RSA procedure such that t σ |σ| is also the arrival time at customer σ i . The algorithm first computes the arrival time t σ |σ| at the last customer when traveling at speed v * fd (minimizing fuel consumption plus driver costs). If this time t σ |σ| is greater or lower than the time windows bounds, it is updated to the nearest bound (Line 6). Then, the necessary speed v ref to arrive at time t σe is computed (Line 7), and the algorithm finds the customer σ p with greatest timewindow violation when using this speed (Lines 10-13). If no violation is found, the solution is returned. Otherwise, the arrival time t σp is updated to the nearest time window bound (Line 15) and SOA is called recursively on two subproblems: from s to p (Line 16), and from p to e (Line 17). Once times t σ i , ∀i = 1, . . . , |σ|, are computed, the associated speeds are revised in such a way that any speed below the optimal speed v * f which minimizes fuel consumption is replaced by a speed v * f and a waiting time (Lines 18-22). The final arrival dates t i at route nodes are obtained (Line 22).
Algorithm 2 is illustrated in Figure 1 on a problem with eight stops. The horizontal axis represents the time, while the vertical axis indicates sequence of customers from bottom to top. The brackets correspond to the time windows and the dots denote the customers' arrival times. Round black dots are related to feasible arrival times; diamond gray dots indicate an early or late arrival time; and square dots indicate the new times after correction. For simplicity, the service times are assumed to be zero. Figure 1a illustrates the arrival times at each customer when the vehicle travels with the optimal speed (v * fd ) starting from σ 1 = 0 at time t σ 1 = a σ 1 . There are two violations and σ 7 Algorithm 2 Speed Optimization Algorithm -SOA if violation > maxV iolation then is the customer with the greatest violation (p = 7). The arrival time at σ 7 is thus adjusted to t σ 7 = b σ 7 (Figure 1b) and two subproblems, P1 and P2, are solved, where P1 considers the subsequence from σ 1 until σ 7 , whereas P2 considers the subsequence from σ 7 until σ |σ| . In P1, the v ref is computed to arrive at σ 7 at time t σ 7 . A new violation is observed at customer σ 4 (p = 4). This violation is corrected ( Figure 1c) and two new subproblems, from σ 1 to σ 4 (P1.1), and from σ 4 to σ 7 (P1.2) are solved. When solving P1.1 and P1.2, no violations are identified (Figure 1c-1d), and thus the recursion is not applied further. This leads to a feasible solution to P1. Subproblem P2 is solved using almost the same rationale of P1 (see Figure1e), except that the customers' arrival times are computed using the optimal speed v * fd (see Line 6 of Algorithm 2). Again, no violation is found, thus completing the whole solution for the problem. At the end, speeds are revised to insert waiting time whenever a speed v < v * f is encountered (Lines 18-22) of Algorithm 2).

Computational Results
In this section we report the computational experiments for the three variants considered in this work, but first we describe the benchmark instances used for evaluating the performance of our algorithm.

Benchmark Instances
The first set of instances considered in this work (Set A) is the one from the PRPLIB, suggested by Demir et al. (2012) and available at http://www.apollo.management.soton.ac. uk/prplib.htm. This set consists of nine different groups, ranging from 10 to 200 customers, each one containing 20 instances. For these instances, the same objective parameter values as However, these instances have a large time windows width, such that it is possible to visit many customers within their respective time windows when traveling at optimal speed as further discussed in Section 5.2. In view of this, we created two additional sets of instances with tighter time windows by modifiying those of the PRPLIB. The time horizon of these new sets is 32400. The time-window width of each customer in Set B is randomly selected between the interval 2000 and 5000 with uniform probability, whereas in Set C it is randomly selected between 2000 and 15000. Once the time windows width are defined, the time window lower bound is randomly selected with uniform probability in the time interval I i which allows to feasibly reach the customer and return to the depot before the end of the time horizon. For a customer i, let W i be the chosen width, and let δ min 0i = d 0i /v max and δ min i0 = d i0 /v max be the minimum driving time from and to the depot, respectively. Then, For the FCVRP and the EMVRP, the well-known benchmark instances of Christofides et al. (1979) and Golden et al. (1998) are used. As in Kara et al. (2007), we use ω = 0.15Q for the EMVRP.

Method performance on the PRP
The proposed matheuristic algorithm was coded in C++ and executed on an Intel Core i7 3.40 GHz with 16 GB of RAM, running under Linux Mint 13. CPLEX 12.4 was used to solve the SP problems. Only a single thread was used, and the algorithm was run 10 times for each instance with different random seeds.
The following parameter values were adopted: n r = 20, n ils = n + 5m, n sp = 150, n pool = 2, T mip = 360 s, and ω tw = 10 8 . The first four parameter values directly follow from  and Subramanian et al. (2013), whereas the latter is a large number to prevent infeasible solutions due to late arrivals. A higher value of T mip than Subramanian et al. (2013) is used, as more time seems to be needed by the SP solver to produce improved solutions.
To investigate the interaction between the local search components, the SP solver, and the speed optimization procedure, we performed further experiments with a static version called ILS-SP-SOA-Stat. For a fair comparison, the CPU time limit for each instance has been set to the average CPU time of the dynamic method for the corresponding instance. In ILS-SP-SOA-Stat, v does not change during an iteration of the algorithm, i.e., the "UpdateSpeedMatrix" function and the Change Speeds perturbation is not used in Algorithm (1). Hence, the task of intelligently recombining routes with different speed distributions is exclusively operated by the SP solver.
The performance of each method is reported in Table 1. Each line corresponds to averaged results on a set of 20 instances. The columns report for each method and instance set the average Gap(%) and CPU time. The Gap(%) for each instance is computed as 100(Z − Z bks )/Z bks , where Z is the objective value of the solution and Z bks is the value of the Best Known Solution (BKS) ever found. Detailed results on each problem instance are provided in B.
It appears that both versions of the proposed algorithm largely outperform the ALNS heuristic of Demir et al. (2012) in terms of solution quality. The new methods also appear to be significantly faster than ALNS, with a speed-up ranging from ×2 to ×50. Note that ALNS experiments have been conducted on a different computing environment and a 3.00GHz CPU. The speed difference between computers of the same generation remains moderate. It can be observed that instances of Sets B and C require more CPU time for convergence, and solution improvements are found more continuously as the search progresses. This is due to tighter time windows. The increased interplay between route and speed optimization leads to a longer sequence of solution improvements, and thus a slower convergence. This aspect is further developed in the next section.
Comparing the results obtained with the two versions of ILS-SP-SOA, the average solutions obtained with the dynamic version are better than those from the static version in all sets of instances but 200-B and 200-C. The dynamic version examines a wider diversity of speed choices. For smaller instances -here up to 100 customers -this allows to converge towards better solutions within the allowed termination criteria. For larger instances with tighter time windows, the benefits of enhanced speed optimization and diversity seems to be counterbalanced by a slower convergence, and the static method, more focused on intensification, may perform marginally better. A longer termination criterion may thus benefit more to the dynamic version that to the static version.

Some insights on instance difficulty for the PRP
To better understand what makes a difficult PRP instance, and how these instances can be better addressed, we further analyze instance characteristics and study how often the optimal speeds are used in the best known solutions. As shown in the following, higher occurrences of "exotic" speeds in the BKS, i.e., speeds that are neither v * f nor v * fd , seem to translate into increased problem difficulty. Figure 2 shows the percentage of the number of arcs in which the vehicle travels at a certain speed in the BKS of each group of instances. Three speed categories are discerned: v * f , v * fd , and other speeds v other . It can be seen that v * fd appears to be used more often than the other speeds in Set A, regardless of the size of the instance. As expected, v * f tends to be used much more often in Set B, which is the one with the tightest time windows. Finally, in Set C, where the time windows are larger than those in Set B, v * fd is the most common in the 10-customer instances, but there is an equilibrium in the 50-, 100-and 200-customer instances between v * f and v * fd . Let %Dist be the percentage of the total distance in which the vehicles are traveling with other speeds than v * f and v * fd in the BKS of a given instance. Figure    A similar analysis is performed in Figure 4 considering the CPU time of ILS-SP-SOA-Dyn, and thus the convergence speed of the method. Overall, the CPU time needed to converge towards a final solution tends to increase for problems with higher %Dist. This suggests that the instance difficultly is correlated with the need of a variety of speeds in the solutions.  This analysis shows that, even if instances with large time-windows are reasonably well solved by methods with separate phases such as the ALNS of Demir et al. (2012) or ILS-SP-SOA-Stat, there is still a notable research effort to be done for more intricate problems, when a larger diversity of speeds must be used to fulfill the time constraints. Because of the variety of antagonist interplays between, speed, distance, and feasibility, more integrated methods that consider jointly speed and route optimization within the neighborhoods are likely to be necessary. ILS-SP-SOA-Dyn is already a first step in this direction, and has demonstrated increased performance on some of the most difficult instances.

Results on other problems
Computational experiments have been conducted on two additional problems with environmental considerations, the FCVRP and EMVRP, in order to assess the performance of ILS-SP-SOA on a wider variety of settings. As previously, a single thread was used and 10 executions were performed for each instance. For a fair comparison with previous methods developed for the FCVRP, we reduced the termination criterion by setting n r = 4, n ils = n/5 + 5m and T mip = 60 s. For the EMVRP we only changed the MIP solver time limit to T mip = 60 s.
For the FCVRP, the results obtained by ILS-SP-SOA were compared with those of Xiao et al. (2012), as can be observed in Table 2. The values highlighted only in boldface indicate that the BKS found by our algorithm equaled the best one reported by Xiao et al. (2012), whereas those in boldface and underline indicate that the BKS found by ILS-SP-SOA improved the one presented by the same authors. To our knowledge, no results on the instances of Christofides et al. (1979) were available for the EMVRP. Hence, Table 3 reports the gap between the average solutions and the BKS found by ILS-SP-SOA.
These results show that ILS-SP-SOA consistently produces solutions of high quality on these two problems. We observe that 22 new BKS solutions were found for the FCVRP, and the average gaps between the average solutions of ILS-SP-SOA and the BKSs is 0.16% for the set of instances of Christofides et al. (1979) and 0.31% for the one of Golden et al. (1998)  do not have a dramatic impact in the solution cost when compared to the associated instances that do not consider these characteristics. Overall, the gaps between the average solutions and the BKS are small, more precisely 0.09% and 0.08%, on the instances of Christofides et al. (1979) and Golden et al. (1998), respectively. When observing such a small variance for VRP problems we can reasonably conjecture that the problem is fairly-well solved.

Concluding remarks
This paper dealt with the Pollution-Routing Problem (PRP), a "green"-oriented variant of the Vehicle Routing Problem (VRP). In order to solve it, we proposed a matheuristic approach, called ILS-SOA-SP, that effectively integrates Iterated Local Search (ILS) with a Set Partitioning (SP) procedure and a Speed Optimization Algorithm (SOA). This approach was also used to solve two other environmental based VRPs, namely the Fuel Consumption Vehicle Routing Problem (FCVRP) and the Energy Minimizing Vehicle Routing Problem (EMVRP). The results of extensive computational experiments demonstrate that ILS-SP-SOA is capable of generating high-quality solutions in a very consistent manner, outperforming previous methods from the literature.
Two new sets of PRP instances were introduced, and our computational experiments show that they are more challenging than the previous one available in the literature. We also showed that the level of difficulty of an instance tends to increase when the vehicles are likely to travel at speeds other than the optimal ones. In other words, good (or even feasible) routing solutions are not easily obtained without using non-optimal speeds, and choosing the correct ones can be a challenging task. The proposed method has the advantage to perform joint route and speed optimization within several local search and integer programming components, and thus performs much better on difficult instances than previously available algorithms.
As a promising line of research, we believe that further efforts are still necessary to incorporate speed decisions more tightly within the local search. This can be done, for example, by developing larger neighborhood structures that are capable of exploring promising regions of the solution space by considering different speeds during the local search. The exploration of these combined neighborhoods is likely to be highly time-consuming, and thus very clever move selection and move evaluation methods may need to be developed.

A Fuel consumption model
The fuel consumption model used in the PRP formulation, given by Eq. (25), is based on the comprehensive emissions model described by Barth et al. (2005), Scora e Barth (2006) and Barth e Boriboonsomsin (2008). It is given by where λ = ξ/κψ and γ = 1/1000η tf η are constants related to the fuel properties, β = 0, 5C d ρA and w are constants associated with the vehicle characteristics and α ij is constant that depends on the road characteristics and on the vehicle acceleration. More precisely, α ij = τ ij + g sin θ ij + gC r cos θ ij , where τ ij is the acceleration and θ ij is the road angle inclination. The meaning and value of all parameters is given in Table 4. Assuming that acceleration and road inclination are null, the fuel consumption can alternatively be written as where w 1 = λkN V , w 2 = λwγgC r , w 3 = λγgC r and w 4 = λβγ. Conversion factor (g/s para l/s) 737 vmin Lower speed limit (m/s) 5.5 (or 20km/h) vmax Upper speed limit (m/s) 25 (or 90km/h) Source:  B Results for the PRP In the tables presented hereafter, Instance denotes the name of the instance, Cost indicates the value of the objective function, |R| is the number of vehicles, CPU T(s) represents the computing time in seconds, Gap(%) is the gap between the solution cost and CPLEX for the 10-customer instances or the best known solution obtained among all experiments for the other instances, given in the column BKS. Note that the approximate optimal solutions for the 10customer instances were computed with the MIP formulation of Demir et al. (2012). It is impractical to consider all possible speed values in the model, and thus a speed-discretization was operated with 500 levels ranging from v * f to v * fd . This leads in general to very high-quality upper bounds.


A Family of Simulation Criteria to Guide DEVS Models Validation Rigorously, Systematically and Semi-Automatically

The most common method to validate a DEVS model against the requirements is to simulate it several times under different conditions, with some simulation tool. The behavior of the model is compared with what the system is supposed to do. The number of different scenarios to simulate is usually infinite, therefore, selecting them becomes a crucial task. This selection, actually, is made following the experience or intuition of an engineer. Here we present a family of criteria to conduct DEVS model simulations in a disciplined way and covering the most significant simulations to increase the confidence on the model. This is achieved by analyzing the mathematical representation of the DEVS model and, thus, part of the validation process can be automatized.

Introduction
The development and use of simulation models has been increasing considerably in recent years. Frequently they are used as the first representation of systems that later will be used for decision-making on critical situations. It has become increasingly necessary the definition of rigorous techniques to assure that these models represent as well as possible the real system being modeled. In other words, Verification and Validation (V&V) of simulation models has become crucial.
According to the U.S. Department of Defense directive [1] verification is "the process of determining that a model implementation accurately represents the developer's conceptual description and specifications". In the context of simulation models, the question that model verification tries to answer is: "Are we simulating, or have we simulated, the model right?". On the other hand, validation is "the process of determining the degree to which a model is an accurate representation of the real-world from the perspective of the intended uses of the model". Here, the question is: "Are we simulating, or have we simulated, the right model?".
Performing V&V of simulation models has been identified as a paramount activity as it can increase the confidence of the user in the simulation results and lead to the accreditation/certification of the simulated system [2]. This accreditation or certification is specially important when the simulation results influence decision making over crucial issues. Figure 1, presented by Robinson [3], which, in turn, is adapted from Sargent [4] shows the different phases of V&V onto the modeling process. Our work is related to the Conceptual Model Validation phase, i.e., validate the conceptual or abstract simulation model against the requirements.

Validation of Simulation Models and Model-Based Testing
The validation of a model against the requirements, usually, cannot be performed mathematically because the requirements are not formal. An alternative way is via simulations. The engineer compares the results of the simulations with the requirements in order to decide if the model is correct or not. This is particularly important when either the model is large, its implementation is critical or critical decisions are made according to the results provided by the implementation of the model or by its simulations. Besides, since it is very important to find as many errors as possible and as earlier as possible, then a thorough simulation process can be an activity that will reduce the total cost of ownership of the target system.
During model simulation it would be desirable to run simulations from all possible simulation configurations and compare these behaviors against the requirements. Unfortunately, exhaustive simulation is impractical in almost all projects, since it involves an infinite number of simulation configurations. Considering this, the selection of an appropriate set of simulation configurations is a crucial issue that should consider two opposite factors: a) the set of simulation configurations should be large enough as to give a reasonable assurance that the model is a correct representation of the requirements; and b) the set should be small enough so V&V fits within time and budget.
In this paper DEVS (Discrete EVent system Specification) [5] is used as the modeling formalism. DEVS has gained popularity in recent years and it is the most general formalism to describe discrete event systems (DES). DEVS is an abstract basis for model specification that is independent of any particular simulation implementation. It is a formalism based on system theory, expressive enough to represent all other DES formalisms, i.e. all models representable in those formalisms can be represented in DEVS [6]. Figure 2 presents a possible DEVS model validation process via simulations. According to this picture, first of all, the requirements and expected results are extracted from the Real World, or from the depiction of the problem being modeled. Based on them, the conceptual or abstract model is defined using a modeling formalism, in this case, DEVS. At this point, it is important to point out that, in order to simplify the validation process, this model must be described using a mathematical or logical representation of DEVS.
According to this validation process, once the abstract DEVS model is defined, a set of simulations configuration is derived from it. Afterwards, both, the simulations and the abstract DEVS model are refined, i.e written in the input language of some concrete simulator (for instance, DEVS-C++ [7], DEVSim++ [8], CD++ [9], PowerDEVS [10], JDEVS [11]). That is, in this context refinement means to translate an abstract representation into a concrete one. The concept of refinement is borrowed from the software engineering community [12] Finally, the (concrete) simulation results need to be abstracted. In this sense, abstraction is the inverse process of refinement. Later, these abstract results are compared with the expected results. This latter issue, generally known as the oracle problem, is acute in this context as it is, for instance, in software testing, and there is no general solution [2,13].
Despite the fact that the final goal is to formalize the whole validation process, this article is focused on the selection of an appropriate set of simulation configurations. As mentioned before, this is the crucial part of this validation process. It consists in turning an infinite problem (the set of all simulations configurations) into a finite one. Further, this must be done trying to keep in the final set the important or revealing simulations. In order words, the final set of simulation configurations must be small and must thoroughly cover all the functional alternatives described in the model. The rest of the validation process consists mainly in an implementation process and is further detailed in Section 4. In general, model simulation is performed according to the experience or intuition of a specialist. Therefore no rigorous guidelines or criteria are followed to define an adequate set of simulations, making the simulation process informal and error prone. On the other hand, given that the selection of the simulations is an informal activity it cannot be automated in a high degree. However, it can be automated to some extent if the selection process is formalized in such a way that, later, a software tool can help in this task. Therefore, it would be desirable to formally define simulation criteria to consider the simulation of a model as a validation process with an acceptable degree of accuracy.
In the software testing field there is an analogous scenario. According to Utting and Legeard [12] software testing deals with the dynamic verification of the behavior of a program on a finite set of test cases, suitably selected from the usually infinite execution domain, against the expected behavior. There are several works trying to formalize the software testing process. Most of these works belong to a subfield of testing known as Model-Based Testing (MBT). Utting and Legeard define MBT as the generation of executable test cases, based on models of the behavior of the system under test (SUT). Hence, an analogy can be established between MBT and model validation by simulation. The model can be seen as the specification of the SUT in MBT and the test case generation as the generation of simulation configurations.
Being so important in software development, the MBT process has been improved up to the point of turning it almost automatic, in many cases obtaining quite good results [13,12,14]. Thus, it is worth to explore if MBT techniques can be replicated in the context of model validation by simulation.
An important part of this automation is possible because there are precise testing criteria. That is, testing criteria indicate which tests must be generated from the model, since there is usually an infinite number of possible tests [12]. Some of these criteria are based on the exploration of the model and others on the exploration of the source code of the program.
Even though many MBT methods can be used or adapted to follow the analogy mentioned before, we based our work on the Test Template Framework (TTF). The TTF is a MBT method presented by Stocks and Carrington [15] and implemented by Cristiá et al. [14]. The TTF was introduced to formally define test data sets providing structure to the testing process. The selection of the TTF as a MBT method is motivated by the fact that it deals with the logical and mathematical definition of the model instead of analyze, for instance, traces or executions of the model. Therefore, it can be applied to systems more complex than Finite State Machines (FSM) and DEVS is much more general than a simple FSM.
Taking all of this into account, the idea of this work is to extend the rigorous validation process of DEVS models presented by Hollmann et al. [16]. Also, a possible way to automate this process is discussed here. Furthermore, a new case study is presented and, in addition, we show some possible errors that could be overlooked during the validation or simulation phase if it is not carried following some systematic and disciplined criteria. The major contribution of this paper is to present an alternative systematic and semi-automatic method to conduct the simulation process of DEVS models in order to validate them. This alternative method is based on well-known techniques in the testing world that can be adapted for the model and simulation community. We believe that simulating DEVS models following the techniques presented in this work will increase the confidence that the model is correct validating aspects or features of the model that could be overlooked by the specialist.
It is important to point out that the issue of running the simulations is not faced in this paper, leaving it as a second step of our research. Here, we give formal criteria to generate the simulation configurations that allow the corresponding simulations to be run in order to validate the model. However, we discuss how this work can be extended in such a way that simulation configurations can be provided to simulation tools.
The remainder of this paper is organized as follows. In the following section we describe and comment some other approaches to model validation. Section 3 presents the simulation criteria that form the core of our contribution. The possible automation of this validation process is addressed in Section 4. In Section 5 two case studies are described, and in section 6 conclusion and some future work is discussed.

Related Works and Similar Approaches
Below we discuss and comment different works involving verification and validation of simulation models. These works either present the most similar approaches to the one presented here, or they motivate or show why our work should be relevant for the modeling and simulation community.
Balci [17] presents guidelines for conducting verification, validation and accreditation of simulation models. He made a classification of the different V&V techniques for simulation models presenting a taxonomy of more than 77 V&V techniques for conventional simulation models and 38 V&V techniques for object-oriented simulation models. In his work, Balci listed model testing as one of the candidates for V&V of simulation models. Also, Labiche and Wainer [2] make a review of the V&V of discrete event system models. They propose to apply or adapt existing software testing techniques to the V&V of DEVS models. In particular, they claim that formal techniques should be applied. Thus justifying our approach. In several works [18,19,20,21] Sargent discusses different approaches, paradigms and techniques related to validation and verification of simulation models. These are interesting works to learn about a generalization of different validation and verification processes, however they do not describe any particular validation process in detail. Our present work complements all these papers by giving a detailed, semi-formal validation method adapted from the software engineering community.
There are several works that use verification techniques, like model checking, to verify the correctness of a model. For instance, Napoli and Parente [22] present a model-checking algorithm for Hierarchical Finite State Machines as an abstract DEVS model. They also focus on the generation of simulation configurations for DEVS, but as counter-examples obtained by the application of their model-checking algorithm. Another relevant and recent work involving verification techniques is [23] where Saadawi and Wainer introduce a new extension to the DEVS formalism, called the Rational Time-Advance DEVS (RTA-DEVS). RTA-DEVS models can be formally checked with standard model-checking algorithm and tools. Further, they introduce a methodology to transform classic DEVS models to RTA-DEVS models, allowing formal verification of classic DEVS. Although model checking techniques are formally defined and they are useful to prove properties and theorems over a model, the main problem of such techniques is the so-called state explosion problem [24], i.e. the exponential blowup of the state space and variables in any real or practical system. This made almost impossible the use of such techniques in large projects, although model checking has been used in real projects.
K. J. Hong and T. G. Kim [25] introduce a method for the verification of discrete event models. They propose a formalism, Time State Reachability Graph (TSRG), to specify modules of a discrete event model and a methodology for the generation of test sequences to test such modules at an I/O level. Later, a graph theoretical analysis of TSGR generates all possible timed I/O sequences from which a test set of timed I/O sequences with 100% coverage can be constructed. Similar to our work, they make an analogy between model verification and software testing.
Another recent work that applies verification techniques over discrete event simulation is [26] where da Silva and de Melo presents a method to perform simulations orderly and verify properties about them using transitions systems. Both, the possible simulation paths and the property to be verified are described as transition systems. The verification is achieved by building a special kind of synchronous product between these two transition systems. They focused their work on the verification of properties by simulation but not on the generation of simulations in order to validate the model.
Li et al. [27] present a framework to test DEVS tools. In their framework they combine black-box and white-box testing approaches. Actually, this work is not really related to ours because they do not validate or verify a DEVS model, whereas they test DEVS implementations. However, it is useful to see how they introduce software testing techniques in the DEVS world.
After reviewing many works about V&V for simulation models it seems that there is no proposals similar to the validation method presented in this article. For instance, we could not find any work that deals with the mathematical or logical representation of the model as a starting point of the validation process. Moreover, we think that, in general, people of the modeling and simulation community do not take this issue into account. Actually they develop their models directly over a simulation tool, using its own modeling language, and make experiments directly over it. By working with the concrete model, these techniques aim at a verification problem, i.e. comparing the concrete model with the abstract model. We, on the other hand, propose, as a complementary technique, to work with the abstract model to check whether it conforms with the user requirements or objectives. More concretely, our work proposes a method to formally validate abstract models and to introduce well known MBT techniques into this community.

A family of simulation criteria
The simulation criteria presented in this work aim to provide a guideline for methodically simulate DEVS models in order to validate them. As we mentioned in the introduction, DEVS is a formalism widely used in the modeling and simulation community and is the most general formalism to describe discrete event systems.
According to Zeigler et al [5], there are two classes of models, Atomic models and Coupled models. Our work concerns only atomic models. In fact, a coupled model is equivalent to a (complex) atomic model [5].
An atomic DEVS Model is defined by the structure: where: • X is the set of input event values, i.e., the set of all possible values that an input event can assume; • Y is the set of output event values; • S is the set of state values; • δ int : S → S is the internal transition function; • δ ext : Q × S → S is the external transition function, where: , s ∈ S, 0 ≤ e ≤ ta(s)} is the total state set, and e is the time elapsed since last transition; • λ : S → Y is the output function; and • ta : S → R + 0,∞ is the time advance function.
δ int , δ ext , λ and ta are the functions that define the system dynamics. Each possible state s ∈ S has associated a time advance, ta(s) ∈ R + 0,∞ , which indicates the time that the system will remain in that state if no input events occur. Once that time has passed, an internal transition is performed, reaching a new state s , s = δ int (s). At the same time, an output event, y, is generated by the output function, y = λ(s). Therefore, δ int , ta and λ define the autonomous behavior of the system.
When an input event arrives, an external transition is performed. The new state depends on the input value, the previous state and also the elapsed time since the last transition. If the system is in the state s and the input event x arrives in the instant e (i.e. e time unites the last transition) the new state, s , is calculated as s = δ ext (s, e, x). In case of an external transition, no output event is generated.

Simulation Configurations Partition
Given a DEVS model, in order to simulate it with some simulation tool, it is necessary to provide an initial state (not defined by the formalism) and a sequence of input events with their corresponding occurrence time, we call this initial state and sequence of input event a Simulation Configuration. Usually, the set of possible simulation configurations (all possible initial state and all possible input events) is infinite, even if the model has finite sets of inputs, states and outputs. Therefore, validation via simulation requires to select some of these possible simulation configurations, see Figure 3. Currently, this selection is done according the experience of the engineer, therefore, a domain expert is usually needed.
The technique presented in this work proposes to divide the set of all possible simulation configurations into equivalence classes by applying one or more partition criteria. We call each of these equivalence classes a Simulation Configuration Class (SCC), i.e. a SCC is a set of possible simulation configurations. Afterward, one simulation configuration of each SCC must be selected, see Figure 4. These selected simulation configurations are the only one that should be executed.
The reason for selecting just one simulation configuration from each class is based on the uniformity hypothesis presented by Bougé et al. [28]. They assert, in software testing, that "A program behaves uniformly on a equivalence class if the following holds: if the program

Simulation Configurations
Partitions Figure 4: Our Proposal: Simulation Configurations Partition works correctly for some input data of the equivalence class then it works correctly for any of them." This is a key assumption made by the software testing community because it allows to reduce the potentially infinite input domain to a small, finite one. Being an assumption it is not proved although many testing methods are based on it [13]. In fact, the very idea of testing is implicitly based on this hypothesis because a single test case actually represents a whole class of test cases. In other words, when a tester selects a test case, he or she is assuming that it represents a set of possible candidates.
We adapt this concept to model validation. Therefore, we say that these subsets, in which the set of possible simulations was divided, are equivalence classes because it is assumed that the model has an uniform behavior for each subset of simulations. Furthermore, if the uniformity hypothesis holds and an error in the model is found for some simulation of a given class, then the same error would be revealed with any other simulation of that class.
In this context, the uniformity hypothesis can be formalized as follows. Let M abs be some abstract model (M abs can be seen as a mathematical or logical formula), M con its concrete model, i.e. its corresponding refinement in a simulation tool language. If a is a simulation configuration derived from M abs , let a be its refinement. Then, M con (a ) means the execution of a on M con ; and M abs (a, M con (a )) asserts that M con (a ) is the expected result with respect to a according to M abs . Note that if M con (a ) is not the expected result of a according to M abs , then M abs (a, M con (a )) if false (i.e. ¬M abs (a, M con (a )) is true). Then, the uniformity hypothesis holds if and only if for every SCC i and a ∈ SCC i the following holds: that is, the uniformity hypothesis holds if the model behaves the same for any two elements of a given SCC i .
With this uniformity hypothesis, some strategies or criteria assume that every element of a class is equivalent to all the others (of the same class) for the simulation of a particular functionality of the model. However, this is just an hypothesis and it is not always satisfied. Therefore, as noted by Stocks and Carrington [15], this assumption is often invalid and they propose to apply repeatedly the strategies in order to partition the classes into sub-classes until either the engineer considers that the classes are reasonable small or each functionality of the model is covered by only one class [29].
Each element of a SCC consists of an initial state (to initialize the simulation) and an input pair containing the event to simulate and its corresponding time, i.e. when the event must be simulated.
Following this idea, in order to describe a SCC we need to define the set of possible initial states and the set of possible input pairs for the simulation. An input pair is an ordered pair (event, time). Therefore a SCC is defined by: • a set of states, IniSt ⊆ S, and • a set of ordered pairs of event and time, where τ represents the no event situation and an internal transitions occur in that case. This is necessary to indicate the simulation of an internal transition.
It is important to point out that in the DEVS formalism no initial state is defined. This, actually, belongs to the simulation phase and, since this work aims to conduct these simulations, an initial state must be defined. This is not necessarily the initial state of the system it is just the starting state for the simulation. Moreover, defining the initial state makes the simulation of the different functionalities of the model simpler, since it would not be necessary to guide the system to a particular state and start a simulation from there.
The total class of simulations, i.e. the set of all possible simulations, for a given DEVS model is defined by: IniSt represents all possible initial states from which a simulation can be started. Note that, initially, all the states of the model, i.e. S, are potential candidates. InP airs is the set of all possible input pairs. Therefore, the idea is to partition IniSt × InP airs by analyzing the DEVS model guided by criteria defined below.

Partition Criteria
Now we present the different partition criteria proposed in this work. The criteria are applied to different aspects of DEVS models. Some criteria apply, for instance, to the external transition function definition, others to the internal transition function definition and others to the definition of the states or input and outputs sets.
It is important to mention that at any time of the partition process it is possible that some criteria could not generate new classes, this is because the result could be classes already obtained. Moreover, it does not matter the order in which the criteria are applied, since each criterion is independent from the others.
In each of the following subsections a partition criterion is described.

Transition Functions Defined by Cases
It is very common to define the external and internal transition functions by cases. The first and more intuitive criterion is to partition the set of possible simulations into several classes, one for each case in the definition of the function.
Let δ ext and δ int be the transition functions of a DEVS model defined by cases: where expr i ext and expr j int are the results of the function if the proposition P i ext or P j int , respectively, holds. This criterion proposes to generate one class for each proposition in the definition of the internal and external transition functions. Each class is defined by: • Associated to the External Transition Function: • Associated to the Internal Transition Function: In the case of the internal transition the idea is to configure a certain state allowing a particular internal transition to occur, that's why the time relative to the τ event is 0.

Extensional Sets
In the definition of DEVS models, sometimes, some sets are defined by extension (states, input events or output events), i.e. listing the elements of the set. Necessarily these sets will be finite and relatively small. Therefore, this criterion proposes to simulate all scenarios where appears at least once each element of these sets.
Let us suppose that the set of states values, S, of some DEVS model is an extensional set: . . , s n } therefore, the classes generated by applying this criterion should be: with i ∈ [1, n]. By applying this criterion the engineer guarantees the simulation of the model in each possible state. However, if this criterion is only applied over the state definition, and no other criterion is applied, certain input events would not be simulated for certain states.
Let us suppose now, that the set of input variables is defined as: the resulting classes now would be: Herein, the engineer ensures to simulate all possible input events. Combining these classes with the former, allows the simulation of all states with all inputs. This is explained further in section 3.3.

Intentional Sets
Set comprehensions or intentional sets is another usual form for defining sets in a DEVS model, i.e. specifying the properties that each element of the set must comply. This is done by a logical predicate, that could be a simple one or a very complex definition involving several operations. Thus, this criterion proposes to use this definition to partition the set of simulation configurations.
Let us suppose that the set of state is an intentional set, S = {s : T Y P E | P (s)}, where T Y P E is the type of the elements of the set and P is a logical predicate. The criterion proposes, first, to write P in its disjunctive normal form (DNF) [30]: . . ∧ P m nm ) and, afterward, partition the simulation configurations set according to this DNF: The same idea can be applied to the set of inputs, outputs or any other intentional set of the model. For instance, let us suppose that the set of inputs X of a given model is X = {(n, m)|n * m > 0 ⇒ n > m}. Therefore, the predicate P = n * m > 0 ⇒ n > m is written is its DNF by applying a known algorithm [30], P = ¬(n * m > 0) ∨ (n > m) and two SCCs should be defined, one for the predicate ¬(n * m > 0) and the other for n > m.

Standard Partitions
In almost all models, different mathematical operators appear in the definitions of the model elements (transition functions, time advance function, state values) and they can be simple (addition, sets union) or more complex (operators defined in terms of simpler ones, functions defined in a programming language or in a pseudo-code). Each operator has a particular input domain and this criterion proposes to divide this domain associating to it a standard partition. A standard partition is a partition of the operator's domain into sets called sub-domains; each sub-domain is defined by the conditions that each operand of the operation must satisfy. Thus, each sub-domain is transformed into a condition to generate a SCC.
Therefore, for each operator in the model a standard partition should be defined. For example, for the operator < (a < b), the standard partition could be [29]: Let us suppose that x 1 , . . . , x n are some of the variables that define the set of states, S, of some model and θ(x 1 , . . . , x n ) is an operator of arity n with the associated standard partition SP 1 (x 1 , . . . , x n ), . . . , SP m (x 1 , . . . , x n ). When θ appears in an expression the set of possible simulations must be partitioned using the standard partition associated with it:

Domain Propagation
This is a particular criterion, since it does not generate new partitions by itself. The purpose of it is to obtain standard partitions of complex operators combining the standard partitions of simpler sub-operators.
Each sub-operation has input domain partitions of its own which are ignored by the standard partitions criterion if it is applied to the complex operator. Using domain propagation the input domain partition of sub-operations are propagated to the higher level [31].
For example, let be a complex operator defined as: (A, B, C) = (A B)♦C where and ♦ are simple operators. Let us suppose that and ♦ have the following standard partitions: We apply first SP to the sub-expression (A B), replacing the formal parameters appearing in SP by A and B respectively: Afterward we do the same with SP ♦ , obtaining: Finally, we combine both propositions obtained and simplify:

Time Partitions
In timed formalisms, like DEVS, modeling the time is a crucial issue. It is very common, in DEVS models, to use additional variables to model time. Furthermore, one characteristic of these models is that the elapsed time appears in the external transition function definition as a variable. Therefore, in the validation of such a model, a question arise: "how do we know if each event has been already simulated in all relevant or significant time interval?". That is, those time intervals in which it is possible to find an error in the model. Again, to answer that question all events must be simulated in all possible time interval making the validation an infinite process.
Therefore, this criterion consists in, first, identifying those variables of the state that interact with the elapsed time or are used to simulate the time advance or timers, for instance. Afterward it defines key time points, or intervals according with the interaction of those variables, for example, using the standard partition for the operations that involve every time variable. Once these key time intervals are defined, one class must be generated for each input at each time point. For example, let us suppose that the time interval [a, b] is relevant to see a particular functionality of the model, therefore it would be meaningful to simulate input events (x, t) with x ∈ X ∪ {τ } and times t < a, t = a, a < t < b, t = b and t > b. Formally, for each defined time interval [a i , b i ], five SCCs should be created: and for each defined time point t j , three classes should be created: The intention of this criterion is to simulate different scenarios where events occur at different moments, to verify the interaction of those variables used to simulate the time among them and the interaction between these ones with the elapsed time (for external transitions).
Later, when this classes are combined with those generated, for instance, by the Extensional Sets criterion, it will be able to simulate all input events in all relevant time instants.

Combining Classes
As has been mentioned in Section 3.2, the partition criteria are independent from each other. Furthermore, each simulation configuration aims to validate a particular characteristic of the model. However, it is usually more efficient to simultaneously validate more than one characteristic. We use efficiency as a measure of the number of errors found in the model. Computational efficiency is not addressed here.
In order to achieve this, it is better to combine the classes generated by applying each criterion. Observe that, however, not all criteria are always applied. This will depend on the model and on the time available for validating it. Also, observe that the same criterion can be applied more than once on the same model. For instance, if the criterion known as Extensional Sets is firstly applied over the state set and secondly over the input set of a given model, the result is two independent sets of configurations. These sets of configurations can be combined in order to simulate the arrival of every input event on every state. In this way, these combined SCCs would find errors, for example, due to the arrival of an event on a wrong state.
Let us see how SCCs can be combined with a toy example. Let M T = (X, Y, S, δ int , δ ext , λ, ta), with X = N and S = N × {ON, OF F }, be part of the model of some system. Suppose that after applying some criteria the following SCCs are obtained: Now, we can combine these classes by intersecting the corresponding IniSt and InP airs sets as follows: InP airs e = InP airs a ∩ InP airs Note, however, that only SCC d and SCC e are valid SCC since SCC f is empty because IniSt f is an empty set. As this example shows, combining SCCs by intersection may yield empty classes. If this is the case: a) eliminate them; and b) keep the SCCs that produced the empty ones.
Observe that, since intersection is commutative, it is the same to combine, for instance, SCC a and SCC b as SCC a ∩ SCC b or as SCC b ∩ SCC a . Moreover, if SCC d is the result of combining SCC a and SCC b , and SCC d is combined with some SCC α the result is SCC b ∩ SCC a ∩ SCC α which is the same regardless of the order in which classes are combined. In summary, SCCs can be combined in any order.

Simulation Sequencing
Once all SCCs are defined it is necessary to set the initial state of the simulation and to execute a transition. However, the new state obtained from the transition may be the initial state of another SCC. If this is the case, then it would be computationally more efficient to continue the simulation by executing the transition indicated by this other SCC. This avoids a new configuration of another simulation for the mentioned SCC, i.e. this reuses the configuration left by the previous run. Then, this process yields a sequence of configurations that should be run one after the other.
Here we propose an algorithm to generate these sequences by analyzing the classes obtained by the application of the criteria. The pseudo code of the algorithm can be seen in Algorithm 1, and it is explained below. select scc ∈ T otSCC 3: select is ∈ scc.IniSt

7:
T otSCC ← T otSCC \ {scc} 8: while ∃ scc ∈ T otSCC | s ∈ scc .IniSt do 9: select scc ∈ T otSCC | s ∈ scc .IniSt 10: select ip ∈ scc .InP airs 11: s ← simulate(s, ip ) SimSeq ← SimSeq ∪ {seq} 16: end while The main progress of Algorithm 1 is led by the following idea: select one SCC (scc ∈ T otSCC), an initial state from that SCC (is ∈ scc.IniSt) and simulate one event of the set of input pairs (event, time) associated to that SCC (ip ∈ scc.InP air). simulate(s, ip ) means to execute on step of the simulation model. Once the selected event has been simulated and a new state, s, is reached the current sequence is updated (sentence 6) and the set of SCCs is reduced by removing the SCC added to the sequence (sentence 7). Now, there are two alternatives: a) wait for an internal transition to occur (if ta(s) = ∞); b) simulate another event. As shown in sentence 8, if there is a SCC, scc , such that s is one of its initial states, then scc is chosen as the next SCC. In which case one of its input pairs is simulated (sentences 10 and 11). Afterward, the process continues repeatedly choosing an event or waiting for an internal transition. At any point of this process, if the state reached from the last simulation step does not belong to any initial state of the remaining SCCs, that sequence ends there and it is added to the set of configuration simulation sequences (sentence 15). A new configuration simulation sequence is started if there are more SCCs.
Thereby, with Algorithm 1 all SCCs are used at least once. Further, more complex coverage criteria could be defined for the generation of the sequences, for example, in a similar way than Souza et al. [32].

Discussion and Automation
As has been mentioned in the introduction, the technique presented in this paper would allow the automation of an important part of the DEVS model validation process. The intention of this section is to show that it would be possible to build a software tool to assist engineers when they apply our method. Therefore, we explain what can be automatically done in each step and we describe the main issues that need to be addressed in order to achieve this automation. Also, further considerations about this methodology are discussed.
An overview of this semi-automated process is the following: 1. Parse the mathematical description of the DEVS model. 2. Select and apply the partition criteria to generate the simulations. 3. Select simulations and generate simulation sequences. 4. Translate the simulation configurations into some simulation language and simulate them. 5. Translate the simulation results into the model formalism. 6. Compare the results with the requirements in order to achieve a verdict about the validation of the model.
The mathematical description of a DEVS model can be automatically be parsed only if it is written in a standard, formal notation. Being so important this issue in the modeling and simulation community, there exists an international group [33] trying to develop standards for a computer processable representation of the mathematical description of DEVS models. This is still an open area. This standard should be based on mathematics and logics, for instance, like Z [34] or TLA [35], instead of a standard based on or akin to a programming language. In this way, engineers do not need to have programming skills for writing their models. For example, we think that DEVSpecL [36], a DEVS specification language developed by Hong and Kim, looks more like a programming language rather than a mathematical representation of DEVS. The definition of a standard mathematical language for DEVS involves the definition of its syntax and semantic and the implementation of a compiler or parser to transform the mathematical model to the input language of a simulator.
Regarding the application of the criteria, a preliminary analysis indicates that it would be appropriate to apply them semi-automatically. We think that a tool should allow engineers to select which criterion must be applied to which part of the model. Moreover, the engineer could add new criteria and use them. Another alternative would be to implement an heuristic that automatically select the criteria according to an analysis over the description of the model.
The application of the criteria involves an important problem, which is the total number of SCCs generated. A preliminary analysis indicates that, although the number of classes could be considerable, it is not exponential since the number of criteria involved is fixed and small. The crucial issue is the combination between classes. The order of classes is given by O(x n ) where x is the number of classes generated by a criterion and n the number of criteria applied (recall that n is fixed and small).
Once the criteria have been applied and the SCCs have been defined in the following phase one simulation configuration for each SCC is selected. Finding a simulation configuration for a given SCC means to find an element belonging to it. Usually, this involves solving a formula over Set Theory, Arithmetic over Integer and Real Numbers and, possibly, other mathematical theories. Further, free variables range over infinite sets, making the problem undecidable. One possibility is to adopt a Satisfiability Modulo Theories (SMT) solver [37,38], by adapting the work of Cristiá and Frydman [39]. Another alternative could be to use constraint solvers such as {log} (pronounced 'setlog') [40,41,42].
The simulation configurations generated above, are described essentially in mathematics. Therefore, to perform the simulations these need to be refined, i.e. rewritten in the input language of some simulation tool, as the model definition. A possible way to do this, is adapting the work done for MBT [43]. This would be a semi-automatic process, since the engineer must define some rules to do this translation.
After performing the simulations, the reverse process must be done. That is, the results of these simulations must be rewritten in the mathematical or formal language used to describe the DEVS model. Again, this would be a semi-automatic process due to the definition of the translation rules. Then, the simulation results can be compared with the expected result (Requirements). This comparison is necessarily manual since involves the requirements.
Besides, it is possible to use this methodology to test the concrete simulation model. However, note that both validations cannot be done simultaneously for the same model. In effect, if our methodology is used to validate the abstract model it is necessary to assume that the concrete model is a faithful representation of the abstract one. On the other hand, if it is used to validate the concrete model it is necessary to assume that the abstract model is correct with respect to the requirements.
Finally, as can be seen, the simulation criteria are applied over the mathematical definition of the model and are based also on mathematical or logical operations. Therefore, the proposed validation methodology is rigorous and systematic (because it is based on mathematics, logic and formal languages). Moreover this methodology is systematic since the whole set of simulation configurations is systematically partitioned obtaining refined and more expressive simulation configuration classes. In this way, given that covering all possible paths of the model is impossible (because they are infinite) we propose to cover the logical and mathematical structure of the model, thus, covering the significant paths of the model.

Case Studies
In this section we show the application of the criteria in two examples. In each case, first we present the requirements of the system, then we describe the DEVS model and finally the SCCs that follow the application of the criteria.

Elevator
The following requirements correspond to the control system of an elevator. It has a control panel with one button for each floor and two other buttons, one for opening and one for closing the door. Furthermore, it has a switch to interrupt or restart its operation. Every time the elevator reaches a floor it receives a signal (the same signal for all floors). To prevent accidents or malfunction of the elevator, it has two sensors, one to check, before or during the door close, if someone or something is crossing it; and the other one to prevent exceeding the weight limit of the elevator. To complete the functionality of the elevator, the system has three timers, whose purpose is described below.
The requirements of the system are: • The door should not be closed if: -Any sensor is active (someone or something is crossing the door, or the weight limit is exceeded).
-The Stop Switch is on.
• If the elevator is stopped at any floor and someone call it from another floor or someone press the button (on the control panel) of another floor, after T D 1 units of time the door should start closing, and takes T D 2 units of time to close completely. However, if the Open Button is pressed or the door sensor activates the timer is reset.
• If after T A units of time (T A > T D 1 ) since the elevator is called the door is not closed, an alarm is fired. Unlike the former timer, this is not reset even though the door sensor is activated or the open button is pressed. It is reset only when the door is totally closed and the elevator starts moving. If the alarm is fired, it should be turned off when this timer is reset. • This elevator has no memory, therefore, it goes to the first floor indicated since it has been stopped. Figure 5, Figure 6 and Figure 7 represent a possible DEVS model corresponding to the control system of the elevator described before.  ((f, f c, eng, d, (ws, ds), sw, a, (at, dt1, dt2, gf t, ot), nt)) =
The input, can be a number (indicating a floor) or different signals, indicating that the elevator has reached a floor, the weight or door sensor has been activated or deactivated, the opening or closing door button is pressed or the switch is turned on or off.
The output, in turn, is a tuple, where each variable represent an indication, respectively, for the display, the engine, the door and the alarm. ⊥ means "do nothing".
Regarding the internal transition function let us describe some cases in order to understand it. For instance, the case (15) represents when the timer of the alarm finish and the alarm must fire, with its corresponding cases (24) and (25) in the output function case. The case (9), in turn, represents when the open button is pressed and the door must be opened, being the door closing. The corresponding case in the output function is (15).
On the other hand, the external transition function is more intuitive to understand each case of it. For instance, case (1) represents when the elevator is called from some floor, different from the actual, the engine is stopped and there is no other floor called. Meanwhile, cases (8) and (9) stand for the case when the switch is activated being the engine stopped or not respectively.

Generating Simulations
Starting with the set of all possible simulations for this example, we now apply each criteria presented before generating different classes. Later, these classes are conjoined, generating new classes. For each criterion, here, we show only some classes, the rest of them can be observed in Appendix A.
Transition Function Defined by Cases. The first criterion that we apply to this example uses the definition of the external and the internal transition functions generating one class for each case in each function definition.
To describe the SCCs for this example, we will use several times a generic s ∈ S defined as s = (f, f c, eng, d, (ws, ds), sw, a, f (at, dt1, dt2, gf t, ot), nt).
Some classes generated by this criterion: Extensional Sets. In this example we have six sets defined by extension, Engine, Door, Sensors, Alarm, Switch and N extT imer. Furthermore, X is the union of an infinite set and a set defined by extension. Since these sets have a relative small number of elements (considering only the finite set of the union resulting in X), we should define one SCC for each element of them, as this criterion proposes: Standard Partitions. Now we apply the standard partitions criterion over the operators <, >, + and −. For these operators could be used the same standard before in subsection 3.2.4. However, some cases must be ignored since the involved variables can not be less than zero (f ≥ 0 and f c ≥ 0). The following classes are some of the result of applying this criterion. The first two relate to the occurrence of the operators < and > in the definition of the internal transition function and the last two, to the occurrence of the operators + and − in the external transition function.
Finally, we apply this particular criterion, taking into account the relation between the elapsed time, e, and the variables used for the timers: at, dt1, dt2, gf t and ot. Considering the values that this variables assume, the relevant key time intervals are: [0, . Therefore, the relevant times t to simulate the each input event are: Hence, some of the classes generated by this criterion are: Combining Classes. Once all criteria are applied we make the intersection of the resulting classes obtaining new and refined ones. Here, we show only some of the classes yielded by these combinations. For example, if we want to test when someone call the elevator from some floor when the engine is stopped with the door open we must conjoin SCC 1 and SCC 49 resulting as follows: • SCC 1 ∩ SCC 49 : InP airs 89 = InP airs 1 ∩ InP airs 49 = Here it can be seen the effect of the commutativity of intersection, since SCC 1 ∩ SCC 49 than SCC 49 ∩ SCC 1 Another interesting combination results from SCC 1 and SCC 57 , because now, the alarm is fired: • SCC 1 ∩ SCC 57 : With the following two classes, errors could be found if the tie-breaking rules are not well defined or they are not given at all. These are obtained by combining SCC 36 with SCC 87 , on one side, and SCC 13 , SCC 59 and SCC 85 on the other. For instance, the SCC defined by IniSt 91 and InP airs 91 simulates the case when the elevator is called at the same time when the "ground floor timer" finish.

Soda can vending machine
This system consists in the control of a soda can vending machine. The machine accepts coins of $ 0.25, $ 0.50 and $ 1. It gives change, optimizing it (i.e. giving the less coins as possible).
The machine has cans of two different prices (normal and diet), and the system that controls the machine must comply with the following requirements: • During an operation, if after T ret units of time no coin is introduced into the machine or no soda is selected, the machine returns all the money that has been introduced.
• Prices of sodas increase as time passes. Every T incr units of time both prices are increased in $ 0.25.
• if the returned money is not collected by the user after T chg units of time the machine recovers it.
• The machine has a display that shows the amount of money introduced or the change after an operation.
• At any time, before selecting a soda, the user can cancel the operation and the machine returns the money.

DEVS Model
In Figure 8 is described a possible DEVS model for this example.
In this model, a state s ∈ S is a tuple where each variable represents, respectively, the machine state, the display, an internal timer, the actual prices (normal and diet), the timer controlling the prices increment, the money stored in the machine, the money inserted for the current operation and the change.
The input values represent, each coin denomination, the request of a normal or diet soda, the cancellation of the current operation and the signal of the change retreated.
The external transition function has one case for each input different from a coin (2, 3, 4 and 5) and one case (1) for all coins. Meanwhile, the internal transition function has one case for each internal state of the machine (m ∈ M achState) for the "operational timer" (1, 2, 3, 4 and 5) and one case (6) for the timeout of the "increase price timer".
The output consists of an ordered pair indicating what to show in the display and how much money (if any) must be returned.

Generating Simulations
As in the previous example, we start with the set of all possible simulations, apply the criteria generating different classes and later conjoin this classes obtaining new ones. As in the previous example, here we show only a few classes. The whole description of the classes obtained can be found in the Appendix B.
Transition Function Defined by Cases. Here we also use a generic s ∈ S to describe the SCCs for this example, defined as s = (m, d, ot, np, dp, it, ms, om, mr).
Some of the classes generated applying this criterion are: • ≥ appears twice (d ≥ np and d ≥ dp) and the standard partition for this operator is equal to the standard partition for the < described in section 3.2.4. The following classes are some the result of applying this criterion. The first two correspond to the application of the criterion over the operation d ≥ np and the last one over the operation d ≥ dp: -IniSt 12 = {s : s ∈ S | d = np = 0}, InP airs 12 = {(x, t) | x ∈ X ∪ {τ }, t ∈ R + 0 } -IniSt 14 = {s : s ∈ S | d = 0 ∧ np > 0}, InP airs 14 = {(x, t) | x ∈ X ∪ {τ }, t ∈ R + 0 } -IniSt 19 = {s : s ∈ S | d > 0 ∧ dp > 0}, InP airs 19 = {(x, t) | x ∈ X ∪ {τ }, t ∈ R+ 0 } • ⊕, by definition, is based on + and with the same type involved, N 0 , therefore they could have the same standard partition. However, since they involve only elements in N 0 no further significant partitions can be proposed. Except if we want to simulate those cases where the implementation of those operations in the modeling language could rise some errors, e.g. overflow errors. In this case, the errors are not properly in the model but in its implementation. This is more related to a testing problem rather than validating through simulations.

Conclusions and future work
We present a family of criteria to conduct DEVS model simulations in a disciplined way and covering the most significant simulations to increase the confidence on the model. The main advantage of performing the simulations of a model as we propose is that users do not need the experience of a specialist or group of specialists, neither a domain expert, to select the simulations to validate the model. The selection of simulations is the result of following a set of formal rules over the mathematical model without having to know about the domain over which it is modeled. This decreases the possibility of overlooking some simulation configurations which could find errors in the model. Another advantage of this work is the possibility of automating at least part of the validation process of DEVS models. An important open issue that needs to be addressed to enable automation is to develop a formal grammar for a mathematical language to describe DEVS models. The development of this grammar is part of our future research.
It should also be considered the possibility of re-using these techniques to test software derived from a DEVS model. A DEVS model could be used as a suitable form of the specification of a system to be implemented in some programming language. The simulation sequences generated from the application of the partition criteria could be used as test cases to test the implementation. Moreover, there exist simulation tools that generate code automatically. Thereby, if the model is thoroughly validated, the resulting piece of software would be correct.
Future work concerns with the automation of the validation process by simulation. First it is necessary to define a standard language to write the mathematical description of a DEVS model. Afterward, we would design and develop the validation tool, including a parser for the standard language, the simulation generator and an automatic translator to some simulation tool. This would also include the definition of new coverage criteria for sequencing simulations.
Other lines of future research are: a) to extend the partition criteria to coupled models; and b) to adapt this validation technique to other formalisms.