4.1. Determining Optimal Credibility level

In order to obtain the optimal amount of Cr parameter the following algorithm which is presented in references [21,22] is utilized.

• Step one: for Cr=0:0.1:1 repeat steps two to five.

• Step two: specify the customers’ demands in the selected Cr level and find the optimal solution to the problem.

• Step three: run the sub program for generating random numbers to clarify the demands of affected areas.

• Sub Program of Generating Random Number

• Sub step one: i=1

• Sub step two: generate a random number of m between the high and low limit of the fuzzy number of the demand of i point and specify the amount of its membership degree (f) based on the membership function of the demand of the affected area i.

• Sub step three: produce the random number of r in [0, 1] interval.

• Sub step four: if r ≤ f, introduce m_i as the real amount of the demand of I point and go to step 4; otherwise, go to step 2.

• Sub step five: increase i one unit. If i ≤ n go to sub step 2; otherwise, go to sub step 6.

• Sub step six: end.

• Step four: with regard to the simulated demand of affected areas calculate the cost of optimized tours at specified Cr level. Consider the cost of occurred failures in calculating the length of tours.

• Step five: do steps three and four 1000 times and calculate the mean of costs.

• Step six: Report Cr content with lowest mean of objective function as the optimized Cr.

5.1. Mathematical modeling

In this section, an integer linear programming model is introduced to formulate the defined problem.

The following assumptions have been considered in this model:

• •

  It is a single-commodity model.

• •

  Each relief center receives goods only from one vehicle.

• •

  Central depot covers no point. In fact, the depot is too far from the affected areas.

• •

  The affected areas receive their needed goods only from one relief center.

• •

  The whole demands of all affected areas must be met.

• •

  The time to travel a distance is the same for all the vehicles and all of them have equal capacity.

• •

  There are some middle points that are not demand points but their distance from the demand points is very little and it is possible to establish temporary depot there.

• •

  It is not possible to establish temporary depot in some affected areas.

The sets considered in this model are as the following:

R′

is the set of all points R′ = {1,2,…R,R + 1}, the point of one is depot and the point of R+1 is virtual.

n

is sub set of R ′ which consists of affected areas (customers of first aid product)

Ba

is the set of affected areas (points with demand) where it is impossible to establish relief centers (barrier areas).

K′

Is the set of vehicles K′ = {1,2,…, K}.

The model parameters are presented below.

t_ij:

The time to travel between two points i and j

d˜i :

Fuzzy demand of point i

dis_max:

Maximum coverage distance

a_ij:

Binary parameter which is equal to 1 if the distance between point i and point j is less than dis_max

Q:

Vehicle capacity

M:

a large number

The variables in the model are as the following:

x_ijk:

is equal to 1 if the distance between two nodes i and j is traveled by vehicle k; otherwise, it is equal to 0

CD_i:

is the collected demands at point i i ∈ {2,3,.., R}

y_i:

is equal to 1 if point i is visited by one of the vehicles; otherwise, it is equal to 0 i ∈ {1,2,.., R}

D_ij:

is equal to 1 if the demand of point i is met by the visited point j; otherwise, it is equal to 0

u_ik:

The time when vehicle k reaches point i

Z:

the time when the good is delivered to the last relief center, in continue, the mathematical model is described as follow:

(8)Minimize Z=max k {u(R+1)k}

(9)∑k=1K∑i=1Rxijk=yj ∀ j=2,…,R

(10)∑j=2Rx1 jk=1  ∀ k∈K

(11)∑k=1K∑i=1R+1xi1k=0

(12)∑i=2Rxi(R+1)k=1 ∀ k∈K

(13)∑k=1K∑j=1R+1x(R+1)jk=0

(14)∑i=1Rxilk−∑j=2R+1xljk=0∀l=2,…,R,k∈K

(15)∑j=2RaijDij=1∀i∈n

(16)∑i=2RaijDij≤Myj  ∀ j=2,…,R

(17)∑i=2Rd˜i Dij=CDj ∀ j=2,…,R

(18)∑i=2R∑j=2R+1CDixijk≤Q ∀ k∈K

(19)yi=Diii=2,…,R

(20)∑i∈Bayi=0

(21)∑k=1Ku1k=0

(22)uik+tij≤ujk+M(1−xijk) ∀i=1,…,R  j=2,…,R+1  k∈K

(23)xijk,yi,Dij∈{0,1}i,j∈R , k∈K

(24)CDi,uik≥0  i∈R , k∈K

Eq. (8) is the objective function of the model which is considered to minimize the time to reach to the last relief center. Constraint (9) states that vehicles only pass through the points where the temporary depots are established. Constraints (10) and (11) indicate that all the vehicles should leave the depot and they cannot return there. Constraints (12) and (13) state that all the vehicles should enter the virtual node (R+1) and no vehicle should leave the node. Constraint (14) states that if a vehicle enters a node, it should leave there. Constraints (15) and (16) state that every affected point should be allocated to a temporary center. Constrains (17) and (18) are related to the capacity of vehicles. Constraint (19) refers to the point that if there are some demands in the visited point, they should not be met by another point within its coverage distance. Constraint (20) guarantees that temporary health center s are not established in inappropriate places. Constraints (21) and (22) are related to omitting the sub tour and calculating the time to reach the temporary health center s. Constraints (23) and (24) identify the type of decision variables.

5.2. Virtual node

In order to model open tours, a virtual node R+1 has been used initiatively. The distance between this node and other points is 0 and all the vehicles after passing the determined point do not return to the origin but go to this virtual node. It is possible to calculate the time to reach the last traveling point using the virtual node. It should be mentioned that the time of receiving to the last point in each tour is equal to the time to reach the virtual node of R+1.

5.3. Model linearization

In the proposed model, the expression (16) is non-linear. In order to linearize expression (8) the expressions (25) and (26) replace the expression (11).

6.1. Harmony Search Algorithm

Due to its applicability for continuous and discrete optimization problems, low mathematical calculations, simple concepts, few parameters and easy implementation, harmony search (HS) algorithm has become one of the appropriate applicable optimization algorithms in various problems in recent years [23]. Compared to the other meta-heuristic methods, this algorithm has fewer mathematical requirements and can be matched with changes in parameters and performances in different engineering problems[24]. In summary, the steps of HS can be stated as the following:

• Step one: Setting the parameters of the problem and algorithm

  In this algorithm, like the other meta-heuristic algorithms, it is necessary to initialize problem and algorithm parameters. HS parameters include harmony memory size (HMS), harmony memory consideration rate (HMCR), pitch adjustment rate (PAR), band width (BW) and number of improvisations (NI).

• Step two: Creating random initial values for harmony memory

  Each harmony vector is randomly initialized and placed in harmony memory. Then, the objective function for each vector is calculated and stored in HM matrix.

• Step three: Creating a new harmony

  In order to produce each component of New Harmony, two different approaches are used in this algorithm. In the first approach which is selected with the same probability as the harmony memory consideration rate (HMCR) a component is randomly extracted from among the corresponding components of one of the harmony memory members. After selecting a component with equal probability to the pitch adjustment rate (PAR), the selected component changes randomly. The range of changes is the same as the band width (BW). In the second approach, the component is randomly initialized with the same probability to (1-HMCR).

• Step four: Updating harmony memory

  If the quality of generated solution is better than the worst solution available in the harmony memory the new solution will replace the worst harmony memory solution; otherwise, the new undesirable solution is eliminated.

• Step five: checking stop criterion

  The third and fourth steps of algorithm continue up to specified number of improvisation and the best answer is reported as the result to the algorithm [23].

  for (j=1 to n) do r1=a random number of continuous uniform of [0 1] if (r1< HMCR) then   Xnew(j) = Xa(j) where a ∈(1,2,…,HMS)   r2=a random number of continuous uniform of [0 1]   if (r2< PAR) then     Xnew(j) = Xnew(j) + r 3 × BW where r3 ∈ (0,1)   end if else     Xnew(j) = LBj + r × (UBj -LBj), where r ∈ (0,1)   end if end for

  Algorithm 1.

  Pseudocode of HS‘s creating new harmony

  Despite the suitability of HS, in some cases it converges quickly. Moreover, this algorithm with a greedy nature is looking forward to improving the worst responses available in harmony memory which in turn causes the diversity of responses in harmony memory to decrease and be trapped by local optimum [24].

  Considering the mentioned difficulties for the HS, the required changes for the development of this algorithm and its efficiency are discussed. In part 6.2 changes applied to HS are explained.

6.2. Improvement of basic HS

In order to improve the quality of HS, four main changes are implied in it as follow.

1. I)

   In new solution generation step in basic algorithm, each component of the new harmony vector is randomly selected from their corresponding components in the harmony memory. However, in new method first a harmony is selected randomly from among harmony memory. For each component of the selected harmony with HMCR probability, it is replaced by the components of one of the members of harmony memory. To do so, with the same probability as HMCR, the rate of component is selected by fitness distance ratio (FDR) which was first introduced by Veeramachaneni et al[25], and then it is determined randomly by 1−HMCR2 and the rest remains unchanged. In FDR the rate of component from one of the answers in memory is selected in such a way that the Eq. (27) is maximized for it [25].

   (27)FDR=f(xi)−f(xj)|xjl−xil|

   In this equation, f(x_i) and x_il are respectively the objective function and the new answer component; f(x_j) and x_jl are the objective function and the selected member of harmony memory.

2. II)

   Pitch adjustment ratio is not randomized and is done in such a way that the new answer components move towards the best available response. After selecting a component randomly, with (PAR) probability it move towards the corresponding component within the best answer available in the harmony memory.

   (28)x=f(best)−f(xi)

   (29)x′=f(best)−f(worst)

   (30)x″=(xx′)×3

   In equations (28) to (30), f (best) and f (worst) are respectively the best and the worst answers available in the harmony memory.

   (31)p=12πe−x″22

   (32)BW=0.1×(1−p)×(xbestj−xij)

   In Eq. (32) x_bestj and x_ij are respectively the rates of j component from the best answers available in harmony memory and selected answer. Given the presented trend for the rate of band width (BW) the more appropriate the selected solution is, a shorter step it takes towards the best solution in the harmony memory, but if the selected solution has an inappropriate objective function it takes longer steps towards the best answer available in harmony memory.

3. III)

   In order to improve the New Harmony, local search is used.

   In this step, the innovative local search between two tours which is explained in 6.4.2. is used to improve the applied response and it is applied to the response as long as the response is not improved.

4. IV)

   In order to place new solution in harmony memory if the quality of new solution is lower than the worst answer in harmony memory, it is eliminated; otherwise, it is added to harmony memory; then S index is calculated via Eq. (33) and the answer with highest rate of S is eliminated.

   (33)S=f(xi)(|f(xi−1)−f(xi)|+|f(xi+1)−f(xi)|)1m

   In the S index f (x_i−1), f (x_i) and f (x_i+1) is the cost function of i-1-th, i-th and i+1-th answers in list of answers, sorted non increasingly. The diversity of available solution in the harmony memory increases through this approach. The flowchart of developed HS is displayed in Fig. 2.

Fig. 2.

flowchart of improved HS algorithm

6.2.1. Solution Representation

In order to represent each solution (harmony) a two-part string with the length of 2R is used where R is the total number of points. In the first part, 0 and 1 are placed. 1 means the establishment of temporary depot and the movement of a vehicle through that point and 0 means no vehicle passes there. If this point has some demand the failure to pass through it means to cover by a passing point. It should be noted that if one of the inaccessible points is visited the solution is modified. In the second string, numbers are placed in the range [0 , 1] that indicate allocation of the points to vehicles. Therefore, the range of [0 , 1] is divided into K (number of vehicles) equal parts with regard to the range to which the number allocated to i section of the second part belongs to. An example of a solution for a problem with 6 places and 2 vehicles is displayed in Fig. 3. According to Fig. 3. nodes 1, 2, 4 and 5 are selected for establishment of temporary depots. In the second part of the string, the numbers less than 0.5 are allocated to vehicle 1 and the numbers more than 0.5 are allocated to the second vehicle. In other words, the demands of points 1 and 5 are met by the first vehicle and the demands of points 2 and 4 are met by the second vehicle.

Fig. 3.

an example of the representing the solution

The method of allocating coverage points to passing points is based on the lowest number of coverage so that if a coverage point can supply its demands through several passing points it will chose a point to which the lowest number of coverage points is allocated.

In each solution, determining the vehicles routs between vertical points is done through innovative algorithm based on GENI which is presented in 6.2.2.

Considering such solution representation, it can be claimed that the constraints of the presented model are evaluated. The only constraints that might be exceeded are the constraints of vehicles capacity. In order to increase the search space, the algorithms are not bound to observe the capacity and exceeding the capacity of vehicles is added to the objective function as a penalty.

6.2.2. Rout construction Algorithm based on GENI Algorithm

GENI algorithm or generalized insertion procedure was first introduced by Gander et al. (1992) [26] for travel salesman problem. This algorithm is generated through the combination of two insertion methods and the results indicate that it has good efficiency in the quality and time of solution.

In this algorithm, at first for selected points, three points that create the greatest coverage are selected and the best path between them is selected so initial tour is created. Then, the following insertion procedure is used for the other points.

Suppose that the initial tour is as (v₁,v₂,…,v_i,…,v_j,…,v_k,…,v₁). Also, if node v is the selected node, N_p(v) is defined as a set of P nodes in the selected tour that have the shortest distance from v compared with the other nodes inside the tour (P is considered to be equal to (3) and P_ij represents nodes between v_i and v_j in the rout of selected tour.

• Sub step 0: Consider node v which is not in the tour.

• Sub step 1: select two nodes v_i, v_j from the set of N_p(v) and the node v_k from the set of N_p(v_j+1)\{v_i,v_j}.

• Sub step 3: Remove edges (v_i,v_i+1),(v_j,v_j+1),(v_k,v_k+1).

• Sub step 4: add edges (v_i,v),(v,v_j),(v_i+1,v_k),(v_j+1,v_k+1) to the tour.

• Sub step 5: Repeat sub steps 1 to 4 for all pairs of the nodes of the set N_p(v)and for all the nodes of N_p(v_j+1)\{v_i,v_j} and report the best state.

• Sub step 6: if there is any node which is not in the tour repeat sub steps 1 to 5.

Finally, two local search methods presented in subsection 6.4 are applied for each solution.

Fig. 4. displays an example of insertion using the above method. Each method continues as long as it improves the solution.

Fig. 4.

An example of GENI algorithm insertion

6.2.4. Adjusting HS parameters

In order to use HS, Taguchi method of experiment designing was used. Taguchi[27] improved a family of matrices of fractional factorial experiments so that he could design experiments after many tests in such a way that the number of tests for one problem decreased. Given that the small kind of objective function of the studied problem (solution variable) is better, the related S/N rate is calculated via the Eq. (35).

(35)S/N=−10 log(1N∑i=1NFi2)

Where, N is the number of test iteration and F_i is the problem solution. The purpose of this experiment is to find the value of harmony search algorithm as the input to achieve optimized response (F) [27].

Harmony memory size (HMS), harmony memory consideration size (HMCR), pitch adjustment rate (PAR) band width (BW), number of algorithm iteration (NI) and the power of factor S (m) are the parameters that this experiment has been used to set. For each one of them 3 values have been considered. In fact, the experiment includes 6 factors and 3 levels. The values considered for each parameter are displayed in Table 1.

parameter	level
	1	2	3
HMS	100	200	300
HMCR	0.5	0.6	0.7
PAR	0.2	0.4	0.6
BW	0.1	0.2	0.3
NI	500	700	1000
M	1	2	3

Table 1.

values considered for harmony search algorithm parameters

According to the number of factors and the levels of each factor, Taguchi plan of L₂₇(3⁶) has been used. In this plan, 27 experiments are considered. By implementing this plan, the experiment of the parameters of basic HS and improved HS algorithm for solving the problem is displayed in Table 2.

parameter	value
HMS	200
HMCR	0.5
PAR	0.6
BW	0.2
NI	1000
m	3

Table 2.

Optimum values of each parameter of harmony search algorithm

6.3. Variable neighborhood search algorithm

Variable neighborhood search (VNS) is a meta-heuristic method or a framework for making heuristic algorithm useful for solving combinatorial problems and global optimization. Its main idea is to change neighborhood by means of local search. Algorithm of variable neighborhood search was first introduced in 1997 by Milladenovic and Hanson [28]. Simple implementation and the quality of solutions resulting from VNS quickly made this algorithm a good method of solving optimization problems.

let NE₁ in which l = {1,2,…l_max} is a predetermined neighborhood structure and NE_l(x) is a set of neighbors of X under the structure of NE₁. VNS algorithm has two main phases of shaking and local search. In the shaking phase, algorithm moves from the current solution (SOL) using l neighborhood structure to another neighbor solution (SOL’). In local search phase, the search is done on Solution SOL’ by means of neighborhood search methods to reach local optimization SOL’*. If the obtained local optimization is better than the current solution, it will replace the solution; otherwise, the next neighborhood optimization (NE₁₊₁) is used to continue the search. The search continues as long as l≤l_max. Pseudo-code of variable neighborhood search algorithm is displayed in the following.

SOL=generate initial solution ();
l = 1;
While (l ≤ lmax)
  SOL′=Shaking (s, NEl)
  SOL′∗=Local search (SOL′)
  if f(SOL′∗) < f(SOL)
    SOL = SOL′∗
    l = 1
  else
    l = l + 1
Until stopping condition are met;
Output The best solution

6.4. Local Search methods

Local search methods improve the created solution. In this research, two local search methods applied to single tour have been used.

6.4.1. 2-opt Algorithm [29]

it is implemented optimally and randomly. Suppose that a tour is defined as j ∈ {1,…, i, i+1,…, j, j+ 1,…, n}. In this algorithm by changing two edges (i, i + 1), (j, j + 1) and adding two edges (i, j), (i + 1, j + 1) a neighborhood is created for the current tour. In optimal algorithm, all possible amounts of i and j (i ∈ {1,2,…,n − 2} and j ∈ {i + 2,…,n}) are examined and the best change is selected in terms of target function; however, in random 2-opt algorithm, I and j are selected randomly and the change in the rout is applied. In the proposed algorithm the optimal algorithm has been utilized. One example of 2-opt algorithm is shown in the Fig. 5.

Fig. 5.

An example of 2-OPT

7.1. Generating Random Samples

For small scale problems the coordinates of each one of the points are determined in the two-dimensional space [10×10] and for large scale problems each of the points is determined in the two-dimensional space [50×50]. Moreover, in all problems, central depot is placed in the origin of coordinates. In order to supply the relative demand of each customer, at first the low limit of fuzzy number is generated by making use of constant uniform distribution with bounds 3 and 8. In order to produce middle limit, the random number has been generated in the range of 1 to 7 and is added to the low limit of demand. Finally, to produce high limit of demand, a number has been generated in the range of 1 to 10 and has been added to the middle limit of demand. The capacity of vehicles is considered in such a way that 0.8 of the total capacity of the vehicles is equal to the sum of high limit of customers’ demands. The range of coverage in small matters is not considered to be fixed and constant, and it is considered as parameter with uniform distribution with parameters between 0.10 and 0.4. In all large scale problems the density of coverage matrix is considered to be 0.2 ± 0.1. According to the rate of density of coverage matrix, the coverage of maximum amount of range can be calculated. The number of points with demand is between 80 to 100% of the total points. The number of points with no demand (not affected from disaster) where it is possible to establish relief operation centers is equal to 5% to 10% of total points. The characteristics of test problems generated are presented in Table 6. and Table 7. in the appendix.

7.2. Determining Optimal Level of Credibility

In order to obtain the optimal level of credibility the algorithm presented in sections 4.1 has been used. First, 10 first small-scale problems and 5 first large-scale problems are selected, (Tables of small and large scale test problems presented in appendix) then, according to the different credibility levels the demands of points are evaluated and according to the estimated demand the optimal routs and related objective function is obtained. In the next step, the tours resulting from the optimal solution are considered and the real demands of the points are supplied using the algorithm of section 4.1 and the objective function of the problem is calculated with regard to the occurred failures. Supplying real demands and calculating the cost are repeated 1000 times and the mean of 1000 times of simulation is considered as the total objective function of the problem at the selected credibility level. The optimal credibility level is determined based on the mean of objective function of total problems. The results are displayed in Table 8. in the appendix. It is worth noting that to solve the small size test problems the CPLEX solver is used and for large scale test problems the improved HS is used

As displayed in Fig. 6., the lowest value of the objective function occurs at the level of 0.6; and so the optimal level of credibility is considered to be equal to 0.6.

Fig. 6.

mean of objective functions and their confidence intervals of generated problems (α = 0.95)

7.3. Efficiency of proposed algorithms for Small-Scale problems

In order to evaluate the efficiency of HS and VNS in small scales, 12 problems according to what was explained in part 7.1, are generated. The data of the problems are available in the appendix.

For the problems generated in small scales, the results of algorithms are compared with the results of optimal solution using CPLEX solver. The summarized results of the comparison are displayed in Table 4., . It is worth noting that for all test problems the optimal level of credibility is considered to be equal to 0.6.

In Table 4. the value of planned optimal objective function is shown by B. In this table, T indicates the running time in second. In the above table, GAP is the percentage difference of optimal planned objective function and the best planned objective function divided by the optimal rate of planed objective function. In other words, the rate of GAP is calculated via the Eq. (36).

In Fig. 7. the solving time of the three algorithms are compared. In terms of solving time, basic and improved HS methods have nearly the same solving time, but in terms of solutions quality improved HS has offered better quality of solution, so that the average error percentage is 1.29% for improved HS and 4.71% for basic HS. The solving time of variable neighborhood search algorithms is relatively less than the other two algorithms but due to the lack of sufficient efficiency the average error of this algorithm is 5.99%.

Fig. 7.

solving time of 3 algorithms for 12 problems in small sizes

The diagram of error percentage of HS and VNS is shown in Fig. 8.

Fig. 8.

error percentage of HS, Improved and VNS algorithms

In order to examine the performance of the proposed algorithms against CPLEX solver, 8 medium size problems were generated and solved by CPLEX and presented algorithm. In all these examples, the time limit for solving by each method is considered to 3600 seconds. In Table 5. the results of numerical examples are provided with an average size.

In Table 5. R.GAP is related gap that GAMS offer after 3600 seconds. The results show that in 9 examples with a median size, CPLEX solver is not able to report the optimal solution for less than 3600 seconds. Finally all algorithms used in the fourth example provided a better solution than the solver CPLEX. As respect to this matter that CPLEX cannot have a better performance than the proposed metaheuristics in median size problems, so in large scale problems we just compare our algorithms.

7.4. Efficiency of proposed algorithms for large-scale problems

In order to investigate the ability of basic HS, improved HS and VNS algorithm in large scale, 12 large scale problems have been generated and solved by metaheuristic algorithms. The data of the problems are presented in the appendix. The summary of the results are illustrated in Table 6., . It is worth noting that for all test problems the optimal level of credibility is considered to be equal to 0.6.

The solving time for large-scale problems of basic and improved HS and VNS algorithm is displayed in Fig. 9. As it is predictable, the mean of VNS solving time is less than that of HS algorithms. The average solving time of basic HS and improved HS is 103.04 and 102.68 seconds respectively while this time for VNS algorithm is 53.54 seconds.

Fig. 9.

solving time for HS and VNS algorithm in large-scale problems

In terms of solution quality, improved HS has been able to offer the mean error of about 0.96%. This value for basic HS is about 3%. The highest rate of solution improvement for the improved HS has occurred in problems 9 and 11 and 12. Fig. 10. shows the rate of improvement by improved HS.

Fig. 10.

improvement through improved harmony search algorithm

Given that comparison between the metaheuristics methods, it is necessary to use nonparametric tests the performance of these algorithms [30]. To this aim the Friedman test was used in the software SPPS. Data from small, median and large-scale are analyzed by Friedman test. In the Fig. 11. output of SPSS software was showed to rank the proposed metaheuristic

Fig. 11.

Results of Friedman test

As can be seen harmony search algorithm developed its lowest rating in terms of objectives and worked better than other algorithm. Also check the information entered in the statistical distribution of error is 0.05 times. Since asymp.sig is less than 0.05. Therefore it follows that average method used for solving the problem are introduced to one another distinction.

In terms of solution quality, HS has found better solutions than VNS algorithm in 9 cases out of 15.

Appendix (A.1)

In Table 7. to Table 10., R is the total number of points and n is the number of points with demands, K is the number of vehicles, Q is the capacity of each vehicle.

In Table 10., the value of planned optimal objective function is shown by B and mean objective function obtained from stochastic simulation shown by M. Also the lower bound and upper bound of confidence interval of mean objective function for each Cr level is reported in the Table 10., . It is worth noting that α is set to the 0.95.