2.1. Fuzzy preliminaries

Since the clustered MTJD involves fuzziness, various fuzzy concepts are discussed briefly in this subsection.

2.1.3 Graded mean integration representation

The graded mean integration representation (GMIR) method is a popular defuzzification technique and has many applications, such as defuzzifying the parameters of an economic order quantity model^[20] and those of the shortest path problem^[21]. Let L⁻¹ and R⁻¹ be the inverse functions of the functions L(x) and R(x), respectively. The graded mean h-level value of the LR-type trapezoidal fuzzy number T˜=(m1,m2,α,β)LR is h[L⁻¹(h)+R⁻¹(h)]/2 as in Fig. 1. Then the GMIR of T˜ is determined by^[22]:

Fig. 1.

The graded mean h-level value of T˜=(m1,m2,α,β).

(4)P(T˜)=∫01h[L−1(h)+R−1(h)]/2dh/∫01h dh

The positive triangular fuzzy number à = (a^l, a^c, a^r) can be regarded as an LR-type trapezoidal fuzzy number à = (a^c, a^c, a^c-a^l, a^r-a^c)_LR. Thus, L⁻¹(x) and R⁻¹(x) of à are given in Eq. (5) and Eq. (6) respectively.

(5)L−1(x)=(ac−al)x+al

(6)R−1(x)=(ac−ar)x+ar

Based on Eq. (4)–Eq. (6), the GMIR of à can be acquired by using Eq. (7)

(7)P(A˜)=∫01h[L−1(h)+R−1(h)]/2 dh/∫01h dh=∫01h[(ac−al)h+al+(ac−ar)h+ar]/2 dh/∫01h dh=∫01(ach2−alh2+alh+arh+ach2−arh2)/2 dh/∫01h dh=12(ac3−al3+al2+ar2+ac3−ar3)/12=al+4ac+ar6

2.2. Representation of the customer satisfaction value

For the distribution of durable goods such as books and consumer electronic devices, many VRP models aim to minimize the distribution costs and use the time window to represent the customers’ time constraints. However, short product lives and fast transportation are two characteristics of the supply chain for fresh and perishable foods^[4]. Chen suggested using an S-shaped curve to depict the relationship between the freshness of perishable food products and profit^[24]. Therefore, a Z-shaped curve is adopted here to depict the relationship between customer satisfaction and arrival time as was suggested in Ref. 12.

Fig. 2 gives a graphical method for illustrating this curve. For example, let a_i = 2 and b_i = 6. When the arrival time of the vehicle to customer i is within two hours, the customer will be very satisfied. As the arrival time passes two hours, the customer satisfaction value (csv) starts to decrease. After six hours, the customer will be very dissatisfied with the delivery.

Fig. 2.

Z-Shaped customer satisfaction function.

To describe the Z-shaped curve mathematically, Eq. (1) is used to depict the customer satisfaction value based on the built-in membership function in Matlab^[25].

2.3. Notations and model formulation

The proposed model consists of the following indices, parameters, and decision variables:

• -

  An index k for the set of m vehicles V = {v₁, v₂, …, v_m}.

• -

  Indices h, i, j for a set of n nodes with customers C = {c₁, c₂, …, c_n}.

• -

  C is divided into u mutually exclusive nonempty subsets C = {C₁,C₂, …, C_u}, each subset representing a cluster of customers.

• -

  C⁺ = C ∪ {c₀} (where c₀ represents the distribution centre).

• -

  Perishable food products are broadly categorized into s types according to their temperature demands for storage. Correspondingly, there should be s types of cabins with different temperature values. The index z is used for the set of temperature values E = {e₁, e₂, …, e_s}.

• -

  Customer demand for z types of perishable food products for customer i is diz.

• -

  The load capacity of each cabin is q.

• -

  The vehicle capacity is denoted by the set B. For example, vehicle k can accommodate B_k cabins.

• -

  pkz cabins are used to store the z type of perishable food products in vehicle k.

• -

  The satisfaction value of the customer i is csv_i, which is generated using Eq. (15). The left extreme and right extreme of the i-th customer satisfaction curve are a_i and b_i respectively.

• -

  Given that the travel time between customer i and customer j cannot be predicted precisely^[12,26], it is considered as a triangular fuzzy number t˜ij=(tijl,tijc,tijr).

• -

  The moment when customer i begins to be served by any vehicle is denoted by the time variable t_i (where t₀ represents the time when vehicles depart from the distribution centre).

• -

  Service time at node i is denoted by w_i.

• -

  M is a large positive number.

• -

  If there exists a vehicle that travels from customer i to customer j (i ≠ j), then xijk=1; otherwise xijk=0.

• -

  If customer i is served by vehicle k, then yik=1; otherwise yik=0.

The mathematical formulation for the clustered MTJD can be stated as follows:

The objective function (16) aims to maximize the sum of the customer satisfaction value. Eq. (17) and Eq. (18) ensure that all customers are served only once. Eq. (19) and Eq. (20) indicate that each route starts and ends at the distribution centre and only one trip is allowed. Eq. (21) states that each customer is served by only one vehicle. Eqs. (22)–(24) are used to represent the relationship between xijk and yik. Eq. (25) guarantees that flow conservation is satisfied. Eq. (26) and Eq. (27) guarantee the consistency of time variables. Eq. (28) indicates that each vehicle accommodates all types of cabins. Eq. (29) guarantees that the total capacity of cabins in any vehicle for accommodating the z type of perishable food products is not exceeded. Eq. (30) ensures that the aggregated capacity of refrigerated containers does not exceed the loading capacity of the vehicle. Eq. (31) guarantees that each customer node belongs to only one cluster. Eqs. (32)–(34) give range constraints on the decision variables.

3.1. Coding method

Each firefly represents a possible solution. For any firefly, its position can be denoted as the set {(v_k, c_i, t_i)|i ∈ C, k ∈ V}. This indicates that the elapsed time when vehicle v_k arrives at customer node c_i is t_i. The light intensity of firefly i is denoted by I_i, which is the total customer satisfaction value.

3.2. Fuzzification and defuzzification

Based on Eq. (7), the GMIR of the fuzzy travel time tij˜=(tijl,tijc,tijr) is defined via Eq. (35).

After the algorithm reads the instance data, tijc is generated via the distance between node i and node j divided by the speed. Next, tijl and tijr are generated based on the value of tijc. The relationship among the three numbers can be determined by the decision maker. Finally, the triangular fuzzy number tij˜=(tijl,tijc,tijr) is defuzzified by applying Eq. (35).

3.3. Population initialization

In many meta-heuristic algorithms including some variants of the firefly algorithm, the initial population is generated randomly^[16,27,28]. Thus, pop_size is defined as the population size, and all fireflies F={F₀, F₁, …, F_{pop_size-1}} are generated at random. Although the initial population of reasonably structured solutions lacks sufficient diversity, there is also a possibility that it could generate high-quality near-optimal values^[29]. For example, Refs. 15 and 17 choose items with high value density to generate initial solutions. Therefore, Algorithm 1 is designed to generate the initial firefly population with structured solutions. Additionally, experiments are carried out with the initial fireflies generated both randomly and by using Algorithm 1 for the purpose of comparison.

for i =0: pop_size-1

  Allocate m neighbouring nodes of the distribution centre from different clusters to each vehicle respectively;

  while not all of the customer nodes are visited

    for j=0: m−1

      if r=rand[0,1] ≤ PI&&all constraints are satisfied

        Allocate the nearest unvisited customer node to vehicle j;

      end if

    end for

  end while

Initialize light intensity Ii;

end for i

To better present Algorithm 1, two examples are shown in Fig. 3. As shown in Fig. 3(a), for node 1, if the random number r is less than or equal to P_I, an edge is generated from itself to its nearest neighbor – node 4. Fig. 3(b) shows the opposite situation. For node 1, if the random number r is larger than P_I, it will give up connecting to its nearest neighbour and wait for the next turn. Thus, node 4 will be connected with node 2 and therefore a solution which differs from that in Fig. 3(a) will be generated. From the above-mentioned presentation, it can easily be inferred that the uniform random variable P_I can prevent the generation of multiple duplicate solutions.

Fig. 3.

Examples of population initialization.

3.4. Movement of fireflies

Given that discrete optimization problems cannot be solved by the original firefly algorithm which was designed for coping with continuous optimization problems, several versions of the discrete firefly algorithm have been proposed^[15–17]. In this algorithm, the arrival time series are used to measure the distance between two fireflies. For example, let ti=(t1i,…,tni) be the arrival time series vector of firefly i. In this sequence, the arrival time of a vehicle to customer h is denoted by thi. Thus, the distance between any two fireflies F_i and F_j is determined by

In the basic firefly algorithm, γ denotes the light absorption coefficient. It is crucial to the convergence speed of the algorithm^[13]. Based on Eq. (36), the movement of firefly F_j toward firefly F_i is determined by

Eq. (37) indicates that the length of the movement of a firefly will be the smallest integral value that is greater than or equal to r_ij/γ. The InsertionFunction is similar to that given in Ref. 16. It selects and extracts one random node from a random route. After that, this node is randomly moved beside a node for which the cluster is same as the extracted one. Fig. 4 shows two examples of how the InsertionFunction works. In Fig. 4(a), node 3 is randomly selected to be moved. Node 11, which belongs to same cluster, is randomly selected as the target node. Next, node 3 is moved to the left of node 11 if the random variable p is less than or equal to 0.5; otherwise it is moved to the right of node 11. As shown in Fig. 4(b), if the extracted node and the target node belong to the same route, the method is also viable.

Fig. 4.

Examples of the InsertionFunction.

Instead of using a repair operator, capacity constraint is taken into account to avoid generating infeasible solutions. This function is repeated n times to enable a movement of firefly F_j. If the new position is superior to the original position, then the firefly is updated.

If F_j is the firefly with the largest light intensity, it moves randomly^[30] in the basic firefly algorithm. Tilahun and Ong pointed out that this behaviour may lead to a reduction in the algorithm performance, because the light intensity of the brightest firefly may decrease in certain directions. They suggested randomly choosing a direction in which the brightness of the brightest firefly increases if the firefly moves in that direction. If there exists no such direction, the brightest firefly must stay in its current position^[31]. However, the brightest firefly F_j will cease moving in the discrete firefly algorithm proposed here, because r_jj = 0 and consequently n = 0 by applying Eq. (36) and Eq. (37). Therefore, n = 1 is set in Eq. (38) for the brightest firefly.

3.5. Repair operation

Sometimes, several customer nodes may exist which are at the end of the route with a csv of zero. The approach mentioned above may generate an impractical route for these, because the visit order does not affect the objective value of the solution. Considering the two routes in Fig. 5 for example, it can be seen that the distance of route 1 is clearly longer than that of route 2, but the sum of csv is the same for the two routes. Therefore, route 1 should be modified by reordering the visiting sequence.

Fig. 5.

An example of the repair operation.

Assume the csv of the successor of node h is zero and the index of the route to be repaired is k. To remedy the problem mentioned above, Algorithm 2 is proposed and presented as follows.

Set J={j|CSV(tj)=0,yjk=1};

Set p = h, yjk=0, xijk=0, tj = 0, ∀i ∈ C, ∀j ∈ J;

repeat

  Randomly select an element s from {j|tpj ≤ tpl, ∀l ∈ J, ∃ j ∈ J};

  Connect node p and node s;

  Calculate ts; Set xpsk=1; Set ysk=1; Set p = s; Set J = J −{s};

until J = ∅

3.6. Procedure for the discrete firefly algorithm

As in the basic firefly algorithm, the movement procedure is repeated until the best solution remains unchanged over rep_gen iterations. To sum up, the pseudocode of the proposed discrete firefly algorithm (DFA) can be summarized according to Algorithm 3.

Set the parameters of the discrete firefly algorithm;

for i=0: n

  for j=0: n

    if i≠j

      Fuzzify tijc and defuzzify t˜ij using Eq. (35);

    else

      tij=0;

    end if

  end for

end for

Initialize the firefly populations;

repeat

  for i=0: pop_size-1

    for j =0: pop_size-1

      if i≠j&& Ii >Ij

        Move Fj using Eqs. 36–38;

      end if

    end for

  end for

  Find Fopt;

  Move Fopt using Eq. (38) (n=1
);

until Fopt remains unchanged over rep_gen iterations

for k=0: m-1

  if a node with zero csv is found in route k

    Repair Fopt using Algorithm 2;

  end if

end for

4.1. Data generation

In the experiments, problems are generated to evaluate the performance of the proposed algorithm. In each instance, the values s=5, q=160L, γ =50, t₀=0, rep_gen=2000 and speed=500m/minute are set. All the customer nodes have corresponding (x,y) coordinates, where both x and y are between −25000m and 25000m. Parameters w_i, a_i, and b_i which refer to time are measured in minutes. In the numerical example about joint distribution in Ref. 10, perishable food products are measured in litres and are categorized into five groups according to their temperature ranges^[10]: (1) tuna, (2) ice cream, ice cubes, and frozen dumplings, (3) milk and juice; (4) biscuits and medicines and (5) lunch boxes and hot meals. Therefore, s is set at a value of 5 and diz is measured in litres in this example. According to general consumer behaviour, customer demand for each type of perishable food product is generated with a uniform distribution, where di1~U(0,5), di2~U(0,8), di3~U(0,10), di4~U(0,15) and di5~U(0,10). As a matter of experience, it may be satisfactory for many consumers to receive the food products within about two-and-a-half hours, and a delivery time of over about six hours could result in customer dissatisfaction. Given that customer attitudes differ with regard to delivery time, the variances of the customer satisfaction parameters are at 900, which can take into account the diversity of customer attitudes appropriately. Therefore, the parameters a_i and b_i are generated with normal distributions N(150,900) and N(360,900) respectively. Given that five minutes or so is enough time for serving a consumer, the parameter w_i is generated with normal distribution N(5,1). Travel times are fuzzified as triangular fuzzy numbers with tijl=0.9tij and tijr=1.1tij. The starting point for each vehicle is (0,0). Given that 10–25 fireflies can deal with most applications^[32], pop_size was set to 20.

4.2. Experimental results

The proposed model and algorithm was applied to 22 instances. Each of these instances was run 10 times. The experimental results are shown in Table 1. DFA₁ was implemented with a randomly generated initial population, and DFA₂ was implemented with an initial population generated by Algorithm 1. Since that the mathematical model of the clustered MTJD is first proposed in this paper, there are no pre-existing dedicated algorithms. Given that Algorithm 1 is a type of greedy algorithm which connects the current node with the nearest unvisited node at each stage when P_I is set to 100%, it can be used as the benchmark for evaluating the effectiveness of the solution approaches, because a greedy algorithm is simple to develop and can solve many applications including the travelling salesman problem and the vehicle scheduling problem^[33].

According to the results presented in Table 1, DFA₁ outperforms DFA₂ in terms of the best value for all the instances except for instances 2, 19, 21 and 22. The median values obtained by DFA₂ are better than those generated by DFA₁ for instances with 120 customer nodes, except for instance 18. All the solutions obtained by DFA₂ are better than those obtained by the greedy algorithm, but the worst value of 10 instances and the mean value of three instances obtained by DFA₁ are even worse than those obtained by the greedy algorithm. Therefore, it can be inferred that DFA₁ is suitable for solving small and medium-sized problem instances. It is also worth mentioning that DFA₂ presents stronger robustness than DFA₁ in general, as shown by the standard deviation values. Additionally, for customers with identical numbers and clusters, more vehicles contribute to higher customer satisfaction values, which partially proves the effectiveness of the algorithms developed in this paper.

Fig. 6 (some dots are overlapping) and Fig. 7 present detailed results of several instances, in order to further investigate the effectiveness of the proposed algorithm. Fig. 8 and Table 2 reveal the efficiency of the algorithm over instances 6 and 16. It can be clearly observed that a larger population size consumes more computational time for both versions of the discrete firefly algorithm proposed in this paper.

Fig. 6.

Graphical description for best operation routes over four instances with pop_size=20.

Fig. 7.

Best vehicle load planning results over No.6 instance.

Fig. 8.

Convergence curves of the discrete firefly algorithms over instance 6.