Fuzzy Uncertainty Management in Multi-Shift Single-Vehicle Routing Problem

Article history: Received: 14 August, 2018 Accepted: 19 October, 2018 Online: 01 November, 2018


Introduction
Vehicle routing problem (VRP) consists in determining a set of routes to visit a fixed set of customers, in order to minimize the total path length. Several versions of the VRP exist. If each customer specifies availability time windows, we deal with a VRP with Time Windows (VRPTW). Basic variants of the VRP consider the route planning for a vehicle fleet in a single period (shift). In this case, vehicles should come back to the depot before the shift ends. This problem is derived from a healthcare routing question. The healthcare company regularly dispatch products to medical sites. In this case, if overtime is allowed, performance could be significantly upgraded [1]- [4]. For instance, if a location scheduled in the next shift is on the current return path to the depot, a limited overtime allows vehicle serving it. This can significantly diminish next shift workload.
Investigation on the connection between health status and work hours documented concern about the influence of working long hours on people fitness [5]. Furthermore, working long hours increases the chance of micro-sleep in car drivers [6]. In these conditions, the company could be accused of vehicle collision due to enormous workload planning [7]. In fact, shift scheduling objective should limit both anxiety and harmful consequences on healthiness maintaining the work shift as steady as possible.
In an effort of limiting extra time, uncertainty effect on scheduling has to be restricted. Frequently, in optimization problems, data are supposed to be known with certainty. Nevertheless, in practice this is rare. Most commonly, the real data depend on uncertainty because of their irregular essence. Because the optimization problem solution usually reveals a great inclination to the data disturbances, neglecting the data ambiguity may conduct to suboptimal or infeasible solutions for a real case. Stochastic VRP (SVRP) was introduced in [8]. See [9] and [10] for a complete review of SVRP.
Robust optimization is a significant technique to deal with optimization problems subject to uncertainty [11]. The main motive for considering robust optimization is the non-stochastic nature of uncertain parameters. In such a case, a methodology is required to analyze compromise between performance and robustness.
Fuzzy set theory is a useful approach to handle uncertain information, while the stochastic approach is suitable to manage the stochastic parameters in VRP [12]. Some VRP papers utilize fuzzy sets theory for studying the influence of uncertain factors [13]- [23].
In this work, the VRP in a fuzzy random context is analyzed and the fuzzy random theory is adopted to manage such uncertain data. Fuzzy random variables illustrate a well-formalized notion involving data obtained from a random test in which data are supposed to be fuzzy sets.
We examine a multi-shift VRP with no overtime allowed considering fuzzy driving and job processing time. The question we considered is derived from a routing problem in maintenance activities. A maintenance team performs jobs in different sites using a vehicle for movements. A crew works in shifts and should come back to the depot before the shift ends. The goal is completing the maintenance activities in various places reducing the makespan. Since we deal with maintenance jobs, we ignore customer waiting times. We investigate the influence of the uncertainty of driving and job processing time on makespan and overtime avoidance.
The body of this paper is structured as follows. In Section 2, we report the problem formulation. In Section 3, we propose the Artificial Immune Heuristic to solve the problem. We present in Section 4 a case study analyzed with the proposed approach. We give concluding remarks in Section 5.
2 Problem formulation 2.1 Classical problem Problem notation is reported in Table 1.
Let N = {1, 2, ..., n} be the set of maintenance jobs to be executed at customer places. Parameter d i,j stands for the travel time between the location of jobs i and j. Parameter q i expresses job i processing time. Parameters e i and l i represent job i time window.
Maximum shift duration is L, whereas p is the number of shifts in the scheduling horizon and P = {1, 2, ..., p} is shift set. We generate p + 1 depot copies described by nodes n + 1, . . . , n + p + 1: • node n + 1 is the origin depot of shift 1, • node n + h represents the destination depot of shift h − 1 together with the origin depot of shift h, with h ∈ 2, . . . , p.
• node n + p + 1 stands for the destination depot of last shift p.
The problem can be represented as a directed graph G = (V , E), where V = N ∪ {n + 1, . . . , n + p + 1}. Arcs (n + h, n + h + 1) are used in the graph to illustrate cases in which a vehicle is not exploited during shift h.
Lastly, b h and c h represent the begin and end time of shift h, with b h = (h − 1)L and c h = hL, h ∈ P .
In this work, the next hypotheses are made. Driving times are positive, i.e. d ij > 0, ∀i, j ∈ V . The triangular inequality is valid for d ij , i.e. d ij ≤ d ik + d k,j , ∀i, j, k ∈ V . At least any single job can be executed in a shift because we suppose shift length satisfies L ≥ d 0,i + q i + d i,n+1 , ∀i ∈ N . Time windows are greater than the activity processing time: l i − e i ≥ q i , ∀i ∈ N . Consequently, a linear formulation for VRPTW is considered. Moreover, a feasible distance set is defined and all the jobs can be processed within their time window and shift duration.
In the following, decision variables are presented: • x i,j = 1 if the crew drives from node i to j, and 0 otherwise, • α i arrival time of the crew at node i ∈ V , i n + 1, • δ i departing time of the crew at node i ∈ V , i n + p + 1, • σ h real shift h duration, • y h = 1 if shift h is used by the crew, and 0 otherwise, The problem is expressed as a mixed integer program (MIP) described in (1)-(13): • The objective function (1) minimizes the system makespan.
• Constraints (2) and (3) provide that the depot nodes {n + 1, . . . , n + p + 1} are always visited. Depot node n + h + 1 is visited after n + h. Depot nodes separate different shifts. All job nodes visited between node n + h and n + h + 1 are served during shift h. If the vehicle is not used in shift h, no job node is visited between node n + h and n + h + 1. Such constraints assure that any shift path starts from the origin depot and ends at the destination depot.
• Constraints (4) ensure that each customer is visited once, that is each node i ∈ N is inspected once.
• Constraints (5) force that for all intermediate nodes (first node n + 1 and last node n + p + 1 are excluded) the inflow is equal to the outflow.
• Constraints (6) are the sub-tour elimination constraints and ensure coherence of time variables. Parameter M is an upper bound of δ i + d ij , ∀i, j ∈ V .
• Constraints (7) represent the connection between arrival and departing times.
• Constraints (13) determine whether shift h is exploited or not. In fact, if x n+h,n+h+1 = 1, then no demand node is inspected between depot h and h + 1 (shift h).

Problem with fuzzy uncertainty
We describe the uncertainty on driving times and job processing times with triangular fuzzy numbers, see Fig. 1. We consider the following notation: • fuzzy travel time: • fuzzy job processing time: this reason, decision variables related to arrival and departing time at node i turn into fuzzy variables α i and δ i . Also, actual shift duration variables σ i become fuzzy. Therefore, fuzzy numbers are introduced both in objective function (1) and constraints (6)- (12).
An important issue arises when models with fuzzy parameters are considered: defining the comparison method for objective function values [24]. In order to solve this question, we have to consider the ranking of fuzzy numbers [25,26,27]. The fuzzy ranking method is part of the solution approach in order to solve mathematical programs having coefficients of the objective function and coefficients of the constraints represented by fuzzy numbers.
The Expected Existence Measure (EEM) operator can be applied to a temporal fuzzy number t, see (14) and [28,29].
Briefly, EEM defines the possibility with which a fuzzy event occurred at a certain time. Assigned a temporal value t 0 , the EEM stands for the possibility a temporal fuzzy number t is lower than t 0 . We denote this as Fig. 1, the membership function µ t (t) of the triangular fuzzy numbert and the corre- Various studies solve fuzzy mathematical problems by exploiting fuzzy ranking methods. First, fuzzy mathematical programming problems are converted into classical mathematical problems. Second, conventional techniques are applied. In the following section, we propose an innovative approach to consider the fuzzy variables in the entire solution process. We present the decision-maker with different optimal solutions based on 2-factor comparison: the objective function value and the possibility degree with which constraints are satisfied. www.astesj.com

Figure 1: Triangular Fuzzy Number
In the previous work [30] we published on this question, we dealt with a different objective function: number of exploited shifts and maximum shift duration. In this research, we consider a different problem in which makespan is minimized. Furthermore, the solution algorithm has been improved with an appropriate solution 'affinity' definition in Section 3.2. Moreover, we extended the experimental campaign in order to investigate our approach performance. For such a reason a mathematical background is introduced in sections 4.2, 4.2.2, and 4.2.3.

Solution Method: Artificial Immune Heuristic
The animal immune system is a versatile patternrecognition system that protects against foreign viruses and bacteria. The immune system is able to recognize and kill germs. The immune system cells, named antibodies, are casually diffused in the body. The system reacts to pathogens and improves the process of recognizing and eliminating pathogens by using two principles: clonal selection and affinity maturation. When a pathogen invades the organism, clonal selection creates a number of immune cells that recognize and eliminate the pathogen. While cellular reproduction occurs, the cells experience high rate physical mutations, as well as a selective process. Cells with superior affinity to the invading pathogen spread into memory cells. Artificial Immune Algorithm (AIA) is a metaheuristic based on such system [31,32,33]. This paper intends proposing a fuzzy artificial immune algorithm to find optimal solutions for the aforementioned problem. AIA notation is reported in Table 2.

Solution Encoding
A solution is encoded as a string by using a fixed-length integer code and providing the order in which nodes are reached. Solution Ψ , representing a scheduling problem with n = 5 jobs on p = 3 available shifts, is reported in Fig. 2. Moreover, Fig. 2 shows the corre-sponding graph path, see section 2. The origin for shift 1 is node n + 1 = 6 that is the path starting node. Then, the crew visits nodes 2, 3 and 5. Shift 1 ends at node n + 2 = 7. After, shift 2 begins and crew inspects nodes 1 and 4. Shift 2 ends when node n + 3 = 8 is reached. Finally, node n + p + 1 = 9 is reached (path end) because no job is executed in shift 3.
Considering fuzzy numbers, we have to appropriately evaluate both the solution objective function f and solution feasibility degree γ.
First, considering the solution Ψ and adopting the fuzzy addition operator described in (15), we calculate actual duration of shift 1, σ 1 , and shift 2, σ 2 , see (16)- (17). Because shift 3 is not used by the crew, we have the same arrival time at node 8 and 9: α 9 = α 8 . We calculate fuzzy makespan as α 9 = α 8 = L + σ 2 . Indeed, makespan is the crew arrival time at the depot at the end of shift 2 (node 8), that is the sum of shift 2 begin time L and actual shift 2 duration σ 2 . Since σ 2 is a fuzzy number, the crisp objective function f in (18) can be adopted using the modal value σ B 2 of the fuzzy number σ 2 .
Second, we determine the degree with which the solution Ψ is feasible. So we calculate the possibility the shift duration constraints are enforced: Φ( σ 1 ≤ L) = γ 1 and Φ( σ 2 ≤ L) = γ 2 . Time windows constraints are considered in the same way. Supposing for solution Ψ , fuzzy shift durations σ 1 and σ 2 are those reported in Fig. 3, it is possible that shift 1 actual duration is greater than maximum limit L. Indeed, half of fuzzy number σ 1 stays on the right side of L. The value γ 1 = EEM σ 1 (L) = 0.5 < 1 indicates the possibility that shift 1 ends before time L. Since, γ 2 = EEM σ 2 (L) = 1,  (19) is determined γ = 0.5.
Eventually, for a given solution, we associate a crisp objective function value f and a feasibility degree γ. We perform a Pareto comparison on (f , γ) pairs for determining the solution Pareto set Γ : f should be minimized and γ should be maximized.
For example, we consider an alternative solution Ψ in which all the available 3 shifts are exploited. The alternative solution Ψ is reported in Fig. 4. Unlike the previous solution Ψ , in Ψ Job 5 has moved from shift 1 to shift 3. On one hand, we have objective function f = 2L + σ B 3 , see (18). On the other, for Ψ , we have γ is equal to 1 because all shift duration constraints are definitely enforced (γ i = 1, ∀i = 1, . . . , N ), see σ C h < 1, ∀h = 1, 2, 3 in Fig. 4. Comparing the solution Ψ to the previous one Ψ , we note that f value gets worse (f Ψ > f Ψ ) and γ values is better (γ Ψ > γ Ψ ). Indeed, with the alternative solution, shifts are completed in time but one additional shift is required. Since we adopt a 2-factor comparison, both solution Ψ and solution Ψ are Pareto optimal.

Affinity
Usually, in meta-heuristic approaches as AIA, an affinity level to be maximized is defined in order to take into account both solution objective value and constraint enforcement. Unfeasible solutions are penalized in terms of affinity. Solutions are ranked on the basis of their affinity values and, eventually, the one with the greatest affinity is defined as 'best solution'. In our approaches, two different factors (f and γ) are both considered in determining the affinity. Solutions having γ = 0 are certainly unfeasible.
We define the solution set as S. Solutions in the Pareto set Γ = {(f , γ)|γ > 0} ⊂ S represent the optimal set, and we assign them the same maximum affinity value as in (20). Note thatf is the makespan of the solution having feasibility degree 1. Since f stands for makespan to be minimized, whereas affinity is to be maximized, we putf in the second term as the makespan upper bound pL occupies the first term.
Using the normalized Manhattan distance from Pareto set (22), the affinity for solutions not in the Pareto set is computed (23). That is, considering a solution (f , γ) not in Γ , the minimum Manhattan distance from Pareto set is calculated and is used as a penalty for affinity. Quantity ∆f is used as a weight in Manhattan distance for balancing the effects of makespan and feasibility degree.

Artificial Immune Algorithm
The proposed AIA is described in the following: Step 1 -Initialization: • (a) Parameter setting: fix the initial population popsize, the number of generations ng, the rate pr1 of Rule1, the rate pr2 of Rule2, for each generation the number of clones nc, the mutation rate mr, the number of mutations nm, the number of exchangeable antibodies nea.
Step 2 -Objective function assessment: • (a) calculate the objective function f in (18) and the feasibility degree γ in (19) for each antibody.
• (c) calculate the affinity for each antibody as reported in (20) and (23).
Step 3 -Clonal selection and expansion: • (a) Take nc antibodies, from the population, with the greatest affinity.
• (b) Produce nc copies of the antibodies considered in Step 3a exploiting a binary tournament rule (randomly select two antibodies from nc antibodies and determine the best antibody).
Step 4 -Generating the next population: • (a) Randomly choose nm antibodies from nc clones and use the mutations to create nm extra antibodies. Apply each mutation operator with the probability mr.
• (b) Include the nm extra antibodies to the current generation.
• (c) Exchange nea worst antibodies with new ones produced like those in Step 1b.
• (d) Copy the Pareto optimal antibodies to the next generation.
Step 5 -Conclusion test: • (a) If the stopping criterion is met, return the Pareto optimal antibodies; otherwise, go to Step 2

At
Step 1a, Rule1 is the full random rule: we randomly select one by one the nodes in the set V \ {n + 1, n + p + 1}. While in Rule2 we choose nodes using a probability that is inversely proportional to the distance between the current node and candidate node. Mutation operator randomly selects two solution indexes (Section 3.1) and swaps content. If the mutation makes unfeasible the depot node sequence n + 2, . . . , n + p, mutation is canceled.
Considering the ordinary AIA approach, the innovation in the method presented in this paper differs in the steps 2b, 2c and 4d. Indeed, step 2b and 2c are used to determine the new antibody affinity. Whereas, step 4d preserves the entire Pareto set in the next generation.

Numerical results
We assessed the performance of the approach described in Section 3.3 by generating different experiments. First, considering a standard non-fuzzy problem, we confronted the performance of the MIP reported in (1)- (13) in Section 2, with the AIA, that is Section 3.3 without steps 2b, 2c, and 4d. Second, introducing fuzzy uncertainty, we adopted the AIA reported in Section 3.3 considering innovative steps 2b, 2c, and 4d.
We developed an application on an x86 family 2.5 GHz Intel Core i7 processor having 4GB RAM and SDD, exploiting "C#.Net4" language. In the following, the AIA parameters and their criteria are defined. Population size pop = 200, number of generations ng = 10000, rate for Rule1 and Rule2 pr1 = pr2 = 0.5, number of clones nc = 20, mutation rate mr = 0.75, mutation number per generation nm = 40, exchangeable antibodies number nea = 20.
Three industrial test cases A, B, and C are adopted. In Table 3 the number of demand nodes n and the number of available shifts p is reported. Maximum shift duration is set to L = 480 min. Since industrial records are protected by a non-disclosure agreement, we report only summary data. The considered geographical is the Salento region in southeast Italy. Traveling and processing times are subject to uncertainty and are set considering user-experience values. Job processing times are equally distributed in three groups: small (around 15 min), medium (around 30 min) and large (around 40 min).

Non-fuzzy version results
Focusing on the non-fuzzy version of case studies A, B, and C, we compared theperformance of the MIP model in (1)-(13) with AIA. Crisp values t B are considered in place of the fuzzy versiont. The non-fuzzy MIP model reported in (1)-(13) was solved with software package IBM CPLEX v.12.5 setting a time limit equal to the running time of the AIA application. In Table 4 we described the achieved results. In particular, for each case study and for each approach, we indicated the time (sec) when a new better solution is detected together with the corresponding solution optimal gap (%). In the first case study A, AIA calculates formerly the MIP optimal solution. For case study B, the optimal solution gap of AIA from MIP is 1% but AIA is tenfold faster than MIP. Finally, in the case study C, the optimal gap is 4%.

Fuzzy version results
Considering the fuzzy version of case studies A, B, and C, we adopted the fuzzy AIA to obtain the solution Pareto set. As reported in Section 3, multiple optimal solutions exist because of the 2-factor comparison. Company manager receives the optimal Pareto set and assesses the best solution. In Section 4.2.1, an example is provided. Assuming different AIA parameter values leads to different Pareto sets as shown in Section 4.2.2.
Our main objective consists in determining the AIA parameters to find as many Pareto set solutions. Consequently, we designed an experimental campaign to determine the best combination of AIA parameter values, see Section 4.

Single experiment results
Considering the fuzzy version of case study A, our software produces the Pareto set reported in Fig. 5. The collected solutions s1, s2, . . . s10, belonging to the optimal Pareto set, are reported in Table 5. For such solutions, jobs are always completed in 3 shifts. Feasibility degree for solution s1 is 1, because σ C i < 480, ∀i = 1, 2, 3. Solution s1 makespan is 2 · 480 + 348 = 1308 min. Since solution s1 is very conservative, the manager is unlikely to accept such a high safety margin.
Whereas, solution s2 has f = 1263 min and γ = 0.968. Indeed, in Table 5 we have σ C 1 = 486 > 480 and σ C 2 = 499 > 480; we note that 96.8% of the area of fuzzy number σ 2 = (349; 424; 499) stays on the left side of 480 (see Fig. 1). The manager prefers solution s2 to s1 as assuming a low risk (3.2%) he/she decreases makespan by 45 minutes. Solutions s3 and s4 have similar feasibility degrees but s4 makespan is lower than s3. Whereas, solutions s5, s6, and s7 have similar makespan but s5 feasibility degree is better than others. Feasibility degrees for solutions s8, s9, and s10 are lower than 0.5: see σ B 2 > 480 min in last three rows in Table 5.
Referring to Fig. 5, manager declares that solutions s2, s4, and s5 are the most valuable and selects the best one based on his/her experience.

Two-experiment results
We compare the experiment performed in previous section 4.2.1 ( = 1), with a new experiment ( = 2) in which pop parameter is decreased from 200 to 100. Consequently, considering the experimental campaign set E = {1, 2}, we report the corresponding solution Pareto sets Γ 1 (pop = 200) and Γ 2 (pop = 100) in Fig. 8. Set Γ 1 , in red, is reported in Fig. 5 too. Set Γ 2 contains 7 solutions, in particular: 4 solutions belong to Γ 1 too, 2 solutions are dominant solutions not included in Γ 1 and one solution is dominated by another in Γ 1 . Instead, one solution in Γ 1 is dominated by a solution in Γ 2 . Since, the new experiment ( = 2) produces two additional solutions in Pareto set and removes one, we have |Γ E | = 11. In Fig. 8, the line represents Pareto front. Calculating ρ 1 = 9 and ρ 2 = 6, we can assess that first experiment has a greater impact than second.

Experimental campaign results
Considering the case study described in Section 4, we designed an experimental campaign. We examined 3 parameters ng, pop, and pr1 on 3 different levels:  Table 6a, experiment impact values are reported for case study A. As parameter pop increases, we obtain a larger number of Pareto set solutions. Increasing the parameter ng from 5000 to 10000 succeeds in increasing ρ, the same cannot be said when passing from 10000 to 15000. Parameter pr1 has a positive influence when pr1 = 0.50.
In Table 6b and Table 6c, experiment impact values are reported for case study B and C. The results confirm the findings of case study A. In conclusion, our solution Pareto set approach significantly depends on population parameter pop. Indeed, in each generation, we preserve the 'entire' solution Pareto set, so it is important having a larger pop value than the non-fuzzy approach. Parameter pr1 = pr2 = 0.50 represents a good compromise between generating Rule1 full random solutions and Rule2 solutions (Section 3.3).

Conclusion
This study presents some real insights into the singlevehicle routing problem with multi-shift and fuzzy uncertainty. Our objective consists of minimizing both the system makespan and shift overtime occurrence. In particular, we investigated the effect of uncertainty in driving and job processing time. We provide optimal solutions for the decision-maker considering a 2-factor comparison: the objective function value (makespan) and the degree with which overtime is avoided. Our approach is currently used by the case study company.