2.1. Problem description

This article aims to address the pharmaceutical distribution routing problem considering customer priority and carbon emissions. The problem can be described as follows.

A particular pharmaceutical distribution center provides cold chain pharmaceutical distribution services to n customers, including hospitals and pharmacies. The hospital-level customers have a higher delivery priority compared to the other customers. The types, locations, time windows, and demands for each customer are known in advance. The pharmaceutical cold chain vehicles are homogeneous, each vehicle has a load capacity of Q, departs from the pharmaceutical distribution center and ultimately returns to the distribution center. Vehicles utilize refrigeration units to maintain the cold chain temperature for pharmaceuticals. During unloading process, heat exchange occurs due to the opening of the vehicle doors, resulting in higher power consumption for the refrigeration unit during the unloading process compared to when the vehicle is in motion. The energy required for the vehicle operation and refrigeration unit operation comes from diesel combustion, which leads to the carbon emissions. The optimization objective of the pharmaceutical distribution routing problem considering customer priority and carbon emissions is to find an optimal solution with minimal cost (including distribution costs, carbon emission cost) and maximal customer satisfaction. Meanwhile, the optimal solution must satisfy a series of complex constraints, for example, the maximum load constraint for vehicles, the demand of each customer must be satisfied, and each customer must be serviced by one and only one vehicle. Furthermore, the main assumptions and constrains of this problem can be described as follows.

1. (1) Only deliver cold chain pharmaceuticals exclusively at customer points, without considering loading at these points.

2. (2) The type, location, time window, and demand quantity information for each customer are all known.

3. (3) Only considering one pharmaceutical cold chain logistics distribution center with known location information.

4. (4) There are sufficient and homogeneous refrigerated vehicles available for distribution. Various types of cold chain pharmaceuticals have the same temperature preservation requirements, and distributing different cold chain pharmaceuticals together in the same vehicle does not affect their quality. The vehicles depart from the pharmaceutical distribution center, complete the day’s distribution, and then return to the pharmaceutical distribution center, forming a closed-loop distribution route.

5. (5) The refrigeration unit operates at its rated power during unloading and at 2/3 of the rated power while the vehicle is in motion (Chen et al., Reference Chen, Liao and Yu2021).

6. (6) The vehicle travels at a constant speed known in advance, without considering factors such as traffic congestion, driver behavior, that may affect the speed and path of the vehicle.

7. (7) Each customer will be served and only by one vehicle.

2.2. Symbol description

The relevant symbols and decision variables used in the pharmaceutical distribution routing model considering customer priority and carbon emissions are shown in Table 1. According to Table 1, one can find the detailed description of these symbols and decision variables.

Table 1. Description of the symbols

2.3. Optimization objective

The optimization objective of the model proposed in this article is to minimize the total cost, including fixed cost, refrigeration cost, fuel consumption cost, carbon emission cost, and penalty cost for violating time window. Each type of cost is detailed as follows.

1. (1) Fixed cost

Fixed cost refers to the associated expenses incurred in using delivery vehicles, including driver salaries, vehicle depreciation cost, and so on. Fixed cost is related to the number of delivery vehicles used in the distribution scheme and can be expressed as the following:

(1) $$ {TC}_1=\sum \limits_{j=1}^n\sum \limits_{k=1}^m{ax}_{0 jk}. $$

1. (2) Refrigeration cost

Pharmaceutical cold chain vehicles use refrigeration units to maintain optimal temperatures for preserving cold chain medication. The refrigeration units consume diesel to generate electricity. Thus, the economic cost of the refrigeration process is assessed based on the cost of diesel consumption. The refrigeration cost consists of two components: the refrigeration cost during the transportation phase and the refrigeration cost during the unloading phase (Chen et al., Reference Chen, Liao and Yu2021). To conserve energy, vehicles shut down their refrigeration units after delivering to the last customer on the route and return to the pharmaceutical distribution center. Therefore, the refrigeration cost while the vehicle is in motion can be expressed as follows:

(2) $$ {TC}_{21}=\sum \limits_{i=0}^n\sum \limits_{j=0}^n\sum \limits_{k=1}^m{x}_{ijk}{P}_d\frac{P_1\left({t}_j^A-{t}_i^L\right)}{E}-{x}_{i0k}{P}_d\frac{P_1\left({t}_0^A-{t}_i^L\right)}{E}. $$

Due to the consistent opening of the cargo compartment doors during the unloading process, there is a sudden increase in the temperature difference between the interior and exterior of the cargo compartment, leading to more intense heat exchange and an increase in refrigeration costs. The refrigeration cost during the unloading process can be expressed as follows:

(3) $$ {TC}_{22}=\sum \limits_{i=1}^n\sum \limits_{k=1}^m{y}_{ik}{P}_d\frac{P_2\left({t}_i^L-{t}_i^A\right)}{E}. $$

Therefore, the total refrigeration cost for the entire distribution process can be expressed as the following:

(4) $$ {TC}_2={TC}_{21}+{TC}_{22}. $$

1. (3) Fuel consumption cost

Fuel consumption cost refers to the cost of fuel consumed by the vehicle engine for travel. This cost is typically associated with the vehicle’s travel distance and payload. In this article, only the fuel consumption cost related to the travel distance is considered and is expressed as

(5) $$ {TC}_3=\sum \limits_{i=0}^n\sum \limits_{j=0}^n\sum \limits_{k=1}^m{x}_{ij k}{P}_d\hskip0.35em \mu {d}_{ij}. $$

1. (4) Carbon emission cost

The refrigeration and travel of delivery vehicles both consume diesel and generate carbon emissions. Therefore, indirect carbon emission costs are introduced (Chen et al., Reference Chen, Liao and Yu2021). The indirect carbon emission cost specifically comprises two parts. The indirect carbon emission cost generated by refrigeration in vehicles can be expressed as

(6) $$ {TC}_{41}=\sum \limits_{i=0}^n\sum \limits_{j=0}^n\sum \limits_{k=1}^m{x}_{ijk}{P}_c\omega \frac{P_1\left({t}_j^A-{t}_i^L\right)}{E}-{x}_{i0k}{P}_c\omega \frac{P_1\left({t}_0^A-{t}_i^L\right)}{E}+\sum \limits_{i=1}^n\sum \limits_{k=1}^m{y}_{ik}{P}_c\omega \frac{P_2\left({t}_i^L-{t}_i^A\right)}{E}. $$

The indirect carbon emission cost generated by the vehicle’s travel can be expressed as

(7) $$ {TC}_{42}=\sum \limits_{i=0}^n\sum \limits_{j=0}^n\sum \limits_{k=1}^m{x}_{ij k}{P}_c\omega \mu {d}_{ij}. $$

Therefore, the total carbon emission cost can be expressed as

(8) $$ {TC}_4={TC}_{41}+{TC}_{42}. $$

1. (5) Penalty cost for violating time window

The timely delivery of pharmaceuticals is crucial to the patient’s safety and well-being. The unit time penalty cost for hospital-level customers is significantly higher than that for pharmacy-level customers. Moreover, the consequences of delivering later than the customer’s time window are more severe than arriving earlier. Exceeding the acceptable threshold results in customer refusal to accept the delivery. Let M represent a sufficiently large positive number, indicating that exceeding the customers’ acceptable time window is not allowed. The relationship between arrival time and penalty cost is illustrated in Figure 1. More specifically, the penalty cost for violating time window for each customer can be expressed as

(9) $$ {p}_i=\left\{\begin{array}{c}{\alpha}_{z_i}\left({ET}_i-{t}_i^A\right),0<{t}_i^A<{ET}_i\\ {}0,{ET}_i\le {t}_i^A\le {LT}_i\\ {}{\beta}_{z_i}\left({t}_i^A-{LT}_i\right),{LT}_i<{t}_i^A\le {LT}_i^{\prime}\\ {}M,{t}_i^A>{LT}_i^{\prime}\end{array}\right.. $$

Figure 1. Penalty cost function for violating time window constraints.

According to Figure 1 and Formula (9), if the time of vehicle arriving at c_i (i.e., t_i^A) belongs to the range of [ET_i, LT_i], the corresponding penalty cost is zero. If the delivery time exceeds the customer’s time window, it incurs a penalty cost that increases linearly with the exceeding time. Moreover, the penalty coefficient for delivery beyond the right time window is higher than that for delivery beyond the left time window, that is, $ {\alpha}_{z_i}<{\beta}_{z_i} $ , and the penalty coefficient for hospital-level customers is higher than that for pharmacy-level customers, that is, $ {\alpha}_1>{\alpha}_2,{\beta}_1>{\beta}_2 $ . Specifically, if 0 < t_i^A < ET_i, the corresponding penalty cost is $ {\alpha}_{z_i}\left({ET}_i-{t}_i^A\right) $ . If LT_i < t_i^A < LT_i’, the corresponding penalty cost is $ {\beta}_{z_i}\left({t}_i^A-{LT}_i\right) $ . Furthermore, the penalty cost for exceeding the acceptable time window (i.e., t_i^A > LT_i’) is set to be M.

Therefore, the penalty cost for violating time window for all customers can be expressed as

(10) $$ {TC}_5=\sum \limits_{i=1}^n{p}_i. $$

2.4. Mathematical model

Optimization model for the pharmaceutical distribution routing problem considering customer priority and carbon emissions is constructed as follows:

(11) $$ \min TC={TC}_1+{TC}_2+{TC}_3+{TC}_4+{TC}_5. $$

Subject to:

(12) $$ \sum \limits_{k=1}^m\sum \limits_{j=0}^n{x}_{jik}=1,i\in C, $$

(13) $$ \sum \limits_{k=1}^m{y}_{ik}=1,i\in C, $$

(14) $$ \sum \limits_{i=0}^n{x}_{ijk}-\sum \limits_{i=0}^n{x}_{jik}=0,j\ne 0, $$

(15) $$ \sum \limits_{j=1}^n{x}_{ijk}=\sum \limits_{j=1}^n{x}_{jik}\le 1,i=0, $$

(16) $$ \sum \limits_{k=1}^m{q}_i{y}_{ik}\le Q, $$

(17) $$ \sum \limits_{i=1}^n\sum \limits_{k=1}^m{y}_{ik}=n, $$

(18) $$ {x}_{ijk}\in \left\{0,1\right\},i\in C,j\in C,k\in V, $$

(19) $$ {y}_{ik}\in \left\{0,1\right\},i\in C,i\ne 0,k\in V. $$

Formula (11) indicates that the optimization objective is to minimize the total cost of five parts, including the fixed cost, refrigeration cost, fuel consumption cost, carbon emission cost, and penalty cost for violating time window, that is, TC ₁, TC ₂, TC ₃, TC ₄, and TC ₅, where TC ₁, TC ₂, TC ₃, TC ₄, and TC ₅ are shown in formulas (1), (4), (5), (8), and (10), respectively. Formulas (12) and (13) indicate that each customer is served only once by a delivery vehicle. Formula (14) indicates that each vehicle must leave after providing service for customer j. Formula (15) indicates that start and end nodes for each vehicle are the pharmaceutical distribution center. Formula (16) indicates that the total demand of customers serviced by the same vehicle must not exceed the load capacity of delivery vehicle. Formula (17) indicates that the distribution center provides service for customers whose number is n. Formulas (18) and (19) represents decision variables.

3.1. Encoding and decoding

This article adopts integer encoding, assigning unique numerical identifiers to each customer, ranging from 1 to n. The length of the chromosome is n + m – 1, where n represents the number of customers, and m represents the maximum available cold chain delivery vehicles for pharmaceuticals. The pharmaceutical distribution center employs identifiers beyond n.

In the decoding phase, the chromosome is divided into multiple paths by extracting identifiers beyond the customer nodes. In each path, the customer nodes are arranged in sequential order based on their positions in the chromosome.

Assuming there are 15 customers and 4 pharmaceutical cold chain delivery vehicles, where customers 1–10 are pharmacy-level customers, and customers 11–15 are hospital-level customers. In the chromosome [6, 3, 16, 11, 7, 17, 14, 8, 5, 15, 1, 2, 4, 18, 13, 9, 10, 12], numbers 1–15 represent customers, and numbers 16–18 represent the pharmaceutical distribution centers, which are utilized to separate different delivery routes. Therefore, this chromosome represents the following routes: path 1: 0–6–3–0, path 2: 0–11–7–0, path 3: 0–14–8–5–15–1–2–4–0, path 4: 0–13–9–10–12–0. Identifier 0 represents the pharmaceutical distribution center. Each delivery vehicle departs from the distribution center, sequentially providing delivery services to the customers along the assigned path, and then returns to the distribution center.

3.2. Population initialization

This article employs a method that a customer is randomly chosen as the starting point. Subsequently, all customers are systematically traversed based on a predefined sequence. Paths are then constructed taking into consideration the vehicle capacity and the left time window. Before adding a customer to the path, check whether the current path’s capacity is sufficient to accommodate the new customer. If the conditions are met, the customer is then inserted into the path at an appropriate position based on the left time window. Once the current vehicle’s load capacity exceeds the constraints, the recorded path is stored, and a new path is initiated. This process is repeated until all customers are allocated to paths. The generated solution will be stored as the first individual in the population.

Following this, an insertion operation is executed on the initial solution. Specifically, two positions are randomly selected, and the gene at one position is inserted into another. After three insertion operations, the resulting individual is added into the initial population. This process continues until a sufficient number of individuals are obtained.

This population initialization process contributes to introducing diversity in the early stages of the GA evolution and provides a promising starting point for the algorithm. This approach facilitates a more effective exploration of the solution space, enhancing the potential of the GA to discover high-quality solutions.

3.3. Crossover with adaptive probability

The select operator employs binary tournament selection, where two individuals are randomly chosen, and the one with lower fitness is eliminated, ensuring the retention of the individual with higher fitness. This iterative process continues until a sufficient number of parent individuals are selected. The use of binary tournament selection enhances the possibility of individuals with superior fitness being preferentially chosen, thereby increasing the probability of their genetic inheritance into the subsequent generation.

To better adapt to the characteristics of the problem, an adaptive crossover probability function is introduced, which flexibly adjusts the crossover probability based on the number of iterations and the fitness level of the population. In the initial stages of the algorithm, a higher crossover probability is employed to promote genetic diversity, ensuring a broad exploration of the search space and increasing the chances of discovering the global optimal solution. As evolution progresses, individuals with higher fitness gradually dominate. At this point, reducing the crossover probability allows the algorithm to focus more on improving convergence speed and the quality of solutions.

The adaptive crossover probability function is defined as follows (Cui et al., Reference Cui, Qiu, Cao, Guo, Chen and Gorbachev2023):

(20) $$ {p}_c={p}_1-\alpha \frac{genF_{lar}}{{ max genF}_{max}}, $$

where P_c is the adaptive crossover probability, P₁ is the base crossover probability, α is the crossover probability adjustment factor, maxgen is the maximum number of iterations, gen is the current iteration number, F_max is the maximum fitness value among all individuals in the population, F_lar is the larger fitness value of the two individuals to be crossed.

3.4. Mutation with adaptive probability

The mutation operator randomly selects either the exchange operator or the insertion operator for execution. The exchange operator randomly swaps genes at two positions. The insertion operator randomly selects two positions and inserts the gene from one position into the other.

The adaptive mutation probability function initially sets a low mutation probability in the early stages of the algorithm to accelerate convergence. As the algorithm progresses, the mutation probability gradually increases to facilitate escaping from local optimum. By dynamically adjusting probabilities, GA can achieve a better balance between the breadth and depth of the search.

The definition of the adaptive mutation probability function is as follows (Cui et al., Reference Cui, Qiu, Cao, Guo, Chen and Gorbachev2023):

(21) $$ {p}_m={p}_2+\beta \frac{genF_i}{{ max genF}_{max}}, $$

where P_m is the adaptive mutation probability, P₂ is the base mutation probability, β is the mutation probability adjustment factor, maxgen is the maximum number of iterations, gen is the current iteration number, F_max is the maximum fitness value among all individuals in the population, F_i is the fitness value of the individual to be mutated.

3.5. Local search based on VNS

VNS is an efficient single solution-based metaheuristic, characterized by strong local search capabilities. The basic idea of VNS is to design multiple neighborhood structures and systematically change the employed neighborhood to escape local optima. More specifically, the VNS starts search from an initial solution, then executes the local search in the employed neighborhood which continually changes according to the search history. The above search process continues until the termination criterion is met (Hansen et al., Reference Hansen, Mladenović and Moreno Perez2010). Researchers have employed this method to successfully address numerous combinatorial optimization problems (Peng et al., Reference Peng, Pan, Gao, Li, Das and Zhang2019, Reference Peng, Deng, Zhang, Shen, Song, Mou and Liu2023; Sun et al., Reference Sun, Shen and Vogel-Heuser2023). In this section, a specialized local search is designed based on VNS to enhance the local search capabilities of the GA. Its primary characteristics are the following two aspects:

1. (1) Three types of neighborhood structures are employed simultaneously (i.e., swap operator, insertion operator, and destruction-repair operator). Search is conducted within each neighborhood structure, if a better solution is found, the algorithm switches to the first neighborhood structure. Otherwise, the next neighborhood structure is set as the current one.

2. (2) In the proposed HGA, local search based on VNS is not applied to all individuals. For certain individuals with lower quality, the local search based on VNS is discontinued. This strategy aims to conserve computational time, achieving a balance between algorithm performance and computational efficiency. More specifically, in addition to the early search phase of the HGA, during the remaining search process, if a solution has an objective function value deviate significantly from the historical best solution up to the current iteration, the application of local search based on VNS is discontinued. For the later search phase of the HGA, if a solution has an objective function value slightly worse than the historical best solution up to the current iteration, the application of local search based on VNS is discontinued.

The specific steps of the local search based on VNS are illustrated as follows:

• Step 1: Set the neighborhood size NSize, individual index i = 1, and neighborhood structure index k = 1.

• Step 2: Decode the ith individual into a delivery scheme and calculate the overall cost of the scheme.

• Step 3: Implement an escaping strategy. If the current iteration exceeds 10% of the allowed maximum iterations and the cost of the current individual exceeds three times of the cost of the best solution found so far, implement an escaping strategy for the ith individual (i.e., not execute the local search), set i = i + 1, go to Step 2. If the current iteration exceeds 90% of the allowed maximum iteration, and the cost of the current individual exceeds 1.5 times of the cost of the best solution found so far, implement an escaping strategy for the ith individual, set i = i + 1, go to Step 2. Otherwise, perform Step 4.

• Step 4: Conduct local search to the ith individual. Repeat the following Steps 4.1–4.3 until the ith individual has executed the following local search NSize times.

  • Step 4.1: If k = 1, utilize the swap operator to generate a neighboring solution. If k = 2, utilize the insertion operator to generate a neighboring solution. If k = 3, utilize the destruction-repair operator to generate a neighboring solution.

  • Step 4.2: If the obtained neighboring solution is better than the ith individual, update the ith individual, set k = 1. Otherwise, set k = k + 1, if k > 3, set k = 1.

  • Step 4.3: Go to Step 4.

• Step 5: Check whether all individuals in the population have completed the above search. If completed, output the optimized population; otherwise, set i = i + 1, go to Step 2.

3.6. Recombination

In binary tournament selection, only one individual is retained from duplicated individuals, and the other duplicated individuals are deleted. Therefore, the number of selected individuals is always less than or equal to the original population size. Following genetic and local search operations on these individuals, the number of offspring obtained still does not exceed the size of the original population. Therefore, the recombination operation retains all offspring individuals obtained after local search in the current iteration, while preserving only the individuals with higher fitness values from the original population.

3.7. Framework of the proposed HGA

• Step 1: Initialization. Set the parameters P ₁, P ₂, α, β, and termination condition. Generate an initial population using the population initialization method in Section 3.2.

• Step 2: Evolution. Repeat the following Steps 2.1–2.6 until the termination condition is satisfied.

  • Step 2.1: Calculate the fitness values of individuals and employ binary tournament selection to choose parents.

  • Step 2.2: Adjust the crossover probability according to the method in Section 3.3, and perform crossover operators on the selected population.

  • Step 2.3: Adjust the mutation probability according to the method in Section 3.4, and perform mutation operators on the population after crossover.

  • Step 2.4: Apply the method in Section 3.5 to conduct local search operations based on VNS on the population after mutation.

  • Step 2.5: Recombination generates a new population by retaining all offspring individuals and the well-adapted individuals from the parental generation.

  • Step 2.6: If the current best individual is better than the best solution found so far, then update the best solution found so far.

• Step 3: If the termination condition is satisfied, then stop the procedure and output the optimization results; otherwise, go to Step 2.

4.1. Experimental instance

A pharmaceutical distribution center needs to provide delivery services to 24 customers, with an ample number of refrigerated vehicles departing from the distribution center. After meeting the cold chain drug demands of all customers, the vehicles return to the pharmaceutical distribution center. Number 0 represents the pharmaceutical distribution center c ₀, and number i represents the customer c_i. The remaining relevant information is shown in Table 2, where customer coordinates, demand, time window, and other data are from Guo (Reference Guo2021), while customer-type data are generated randomly.

Table 2. Information of the pharmaceutical distribution center and customers

Moreover, the parameter settings required for the proposed optimization model are shown in Table 3. Furthermore, about the algorithm parameters of the proposed HGA, the population size is set as 50, the maximum number of iterations is 250, the base crossover probability is 0.9, the base mutation probability is 0.05, and both the crossover and mutation probability adjustment factors are set as 0.5.

Table 3. Parameters of the model

4.2. Comparison results with four optimization algorithms

To validate the effectiveness of the presented HGA, we compare it with four optimization algorithms: improved genetic algorithm (IGA; Cao, Reference Cao2021), GA, improved ant colony algorithm (IACA; Guo et al., Reference Guo, Gao, Shi, Wu, Wang and Zhang2023), and simulated annealing (SA). Each algorithm is run independently 20 times on the same test case. The experiments are performed on a computer with Intel(R) Core (TM) i7-7500U CPU @ 2.70 GHz 2.90 GHz, 8.00 GB RAM, and Microsoft Windows 10. All the algorithms are coded in MATLAB 2020b.

To ensure equal comparison among the algorithms, identical parameters are employed. The remaining parameter settings for the IGA and GA are as follows: the population size is set to be 50, the maximum number of iterations is set to be 250, the crossover probability is set to be 0.9, and the mutation probability is set to be 0.05. The remaining parameter settings for the IACA are as follows: the ant quantity is 50, the pheromone heuristic factor is 1, the expected heuristic factor is 5, the pheromone evaporation factor is 0.5, and the maximum iteration number is 250. The remaining parameter settings for the SA are as follows: the initial temperature is 1,000, the cooling rate is 0.99, the maximum iteration number for the inner loop is 50, and the maximum iteration number for the outer loop is 250.

The best and average cost obtained from 20 runs of each algorithm are shown in Table 4, which represent the best value and the mean value achieved in the 20 runs, respectively. The HGA designed in this article consistently exhibits superior performance, with a best cost of 986.92 yuan and an average cost of 990.18 yuan, both surpassing those achieved by other algorithms. The IGA inferior exhibits noteworthy performance, with a best cost of 992.85 yuan, slightly higher than that attained by the HGA. The IACA attains a best cost of 1,268.54 yuan, while the GA attains a best cost of 1,878.48 yuan, and the SA achieves a best cost of 2203.23 yuan. However, its relatively higher average cost reflects shortcomings in terms of stability. In contrast, both the GA and SA exhibit certain limitations in addressing the optimization problem presented in this article.

Table 4. Test result of different algorithms

The best path optimization solution obtained among the 20 runs of the HGA is shown in Table 5, and the corresponding distribution route is displayed in Figure 2. The total cost of this solution is 986.92 yuan, with the fixed cost of 150 yuan, the refrigeration cost of 237.66 yuan, the fuel consumption cost of 561.72 yuan, the carbon emission cost of 37.54 yuan, and no penalty cost for violating the time window. The delivery scheme utilizes three vehicles, following the delivery routes: Route 1: 0–8–22–2–1–11–24–15–18–0, Route 2: 0–14–17–13–4–7–23–12–0, Route 3: 0–21–10–9–6–5–3–16–19–20–0, where the number 0 represents the pharmaceutical distribution center. In conclusion, the proposed HGA is an efficient approach for dealing with the pharmaceutical distribution routing problem considering customer priority and carbon emissions.

Table 5. Path optimization results obtained by the proposed HGA

Figure 2. Best route obtained by the proposed HGA.