Model and algorithm for pharmaceutical distribution routing problem considering customer priority and carbon emissions

Abstract

Pharmaceutical distribution routing problem is a key problem for pharmaceutical enterprises, since efficient schedules can enhance resource utilization and reduce operating costs. Meanwhile, it is a complicated combinatorial optimization problem. Existing research mainly focused on delivery route lengths or distribution costs minimization, while seldom considered customer priority and carbon emissions simultaneously. However, considering the customer priority and carbon emissions simultaneously will not only help to enhance customer satisfaction, but also help to reduce the carbon emissions. In this article, we consider the customer priority and carbon emission minimization simultaneously in the pharmaceutical distribution routing problem, the corresponding problem is named pharmaceutical distribution routing problem considering customer priority and carbon emissions. A corresponding mathematical model is formulated, the objectives of which are minimizing fixed cost, refrigeration cost, fuel consumption cost, carbon emission cost, and penalty cost for violating time windows. Moreover, a hybrid genetic algorithm (HGA) is proposed to solve the problem. The framework of the proposed HGA is genetic algorithm (GA), where an effective local search based on variable neighborhood search (VNS) is specially designed and incorporated to improve the intensification abilities. In the proposed HGA, crossover with adaptive probability and mutation with adaptive probability are utilized to enhance the algorithm performance. Finally, the proposed HGA is compared with four optimization algorithms, and experimental results have demonstrated the effectiveness of the HGA.

Impact Statement

Pharmaceutical distribution routing problem considering customer priority and carbon emissions is an extension of pharmaceutical distribution routing problem by further considering customer priority and carbon emissions simultaneously. The pharmaceutical distribution routing problem considering customer priority and carbon emissions is an important problem, since efficient schedules will not only reduce operating costs and carbon emissions, but also enhance customer satisfaction. To address the problem, this article establishes the corresponding mathematical model, and further devises an effective hybrid genetic algorithm (HGA). The performance of the devised HGA is enhanced by developing and merging an efficient local search based on variable neighborhood search (VNS). Comparison results have shown the proposed HGA is an effective approach for solving the aforementioned problem.

1. Introduction

With the continuous breakthroughs in medical technology, biomedical science has afforded novel opportunities for the treatment of severe and rare diseases. The rising pharmaceutical expenditure has become a pressing concern for many nations (Brekke et al., 2023). In this context, the significance of the pharmaceutical industry is increasingly pronounced. Within the development of pharmaceutical industry, the pharmaceutical distribution routing has become a critical issue. Effective optimization of pharmaceutical distribution ensures timely access to appropriate medical resources for patients, enhances the economic efficiency of healthcare institutions, and improves the overall operational efficiency of the healthcare system. Currently, pharmaceutical logistics confronts challenges such as the continuous expansion of pharmaceutical varieties, stringent transportation requirements, and the escalating societal demands. Hence, addressing pharmaceutical distribution routing problem through scientifically effective path planning has become a pressing issue that requires solution.

The pharmaceutical distribution routing problem is a typical NP-hard problem, posing computational complexity in finding solutions (Nakiboglu and Gunes, 2018). The existing solution approaches primarily encompass mathematical programming methods and metaheuristic techniques. In order to ensure the secure management of health-care waste, Ouertani et al. (2023) proposed a transportation model that segregates hazardous waste from nonhazardous waste. They designed an adaptive genetic algorithm (GA), demonstrating its efficiency through experimental study on benchmark instances and real data instances. A healthcare facilities clustering method applicable to hospitals in Moroccan was introduced by El Midaoui et al. (2021), who proposed a simulation of the multi-depot vehicle routing problem (VRP) solved using a GA. GA was used effectively by Zheng et al. (2011) to address the optimization problem of dynamic drug distribution routes in urban areas. Yan (2020) improved the delivery efficiency and resource utilization for a pharmaceutical enterprise by applying improved GA and saving method to design distribution scheme. Bouziyane et al. (2020) studied the disrupted VRP with soft time windows. They proposed an improved multi-objective local search utilizing a combination of GA and variable neighborhood search (VNS) to improve the efficiency of pharmaceutical distribution. Moghadam (2010) studied the pharmaceutical distribution routing problem under demand uncertainty and unknown distributions, successfully using a particle swarm optimization (PSO) to address real world case. Redi et al. (2020) presented a simulated annealing (SA) heuristic to finding a set of vehicle distribution routes with minimal transportation time for a pharmaceutical company in Jakarta. Pacheco and Laguna (2020) focused on the mask distribution routing problem during the COVID-19 pandemic. Their model aimed to minimize both the time of the longest route and the total travel distance, with a corresponding tabu search (TS) designed for solution. Zhu and Wang (2021) used ant colony optimization (ACO) to solve the distribution model considering traffic conditions and the distribution model without considering traffic conditions. The comparative analysis indicated that considering traffic conditions in pharmaceutical cold chain distribution results in a 11.96% reduction in overall costs.

Petroianu et al. (2021) considered the challenges of insufficient vehicles, road availability, and weather conditions in both rural and urban impoverished areas for vaccine distribution, and developed a path optimization tool named route optimization tool (RoOT) based on index algorithms. This tool was designed to facilitate the efficient distribution of vaccines and other medical supplies to healthcare centers across various regions in Mozambique. Zhao et al. (2021) adopted a space–time network model to describe the decision-making process and simplify the model structure. They designed a comprehensive optimization algorithm based on the Lagrangian relaxation method, effectively addressing the challenging problem of drug collaborative distribution. Campelo et al. (2019) developed an instance size reduction algorithm and a mathematical programming-based decomposition approach to address the consistent VRP and successfully applied to real-life scenarios. Taking into consideration, the actual demand for medical supplies such as masks and vaccines during the COVID-19 pandemic, researchers divided the distribution into two sub-problems. They addressed the distribution challenges in major cities by combining machine learning and optimization algorithms, offering scheduling guidance for similar issues (Chen et al., 2022; Zheng et al., 2023). A hierarchical decision mixed-integer linear program was proposed by Guo et al. (2023), who gave a case study of COVID-19 vaccine administration to demonstrate improved efficiency in vaccine administration and distribution, along with a reduction in overall costs. Inspired by the production and distribution challenges in a traditional Chinese medicine company, Qiu et al. (2023) proposed a mixed integer linear programming model that considers hybrid flow shop production and multi-trip, multi-vehicle distribution. They developed an improved memetic algorithm (MA) to address this joint optimization problem and substantiated the superiority of integrated scheduling through numerical experiments.

The existing pharmaceutical distribution routing problem mainly focuses on planning the shortest or lowest-cost routes, with limited research considering the differentiated demands of customers. Customers in pharmaceutical distribution include hospitals, retail pharmacies, and other entities. Among them, hospitals directly provide medical services to patients, and the demand for medications is typically urgent. Therefore, it is crucial to ensure that hospitals receive high priority service during the distribution process. The research on pharmaceutical distribution routing problem considering customer priority has significant value. With the increasing delivery requirements of customers, considering customer priority has become a novel and crucial factor. It allows for the effective segmentation of customers into different levels, allocating corresponding service levels based on their specific demands and importance. This differentiated delivery service not only enhances customer satisfaction but also contributes to boosting the competitiveness of the business in the market.

In the existing research on VRP, researchers have gradually recognized that customer priority not only affects service quality but is also closely related to resource utilization and economic efficiency. Studies have been conducted to address VRP considering customer priority. Doan et al. (2021) classified customers into different priority groups to ensure that high priority customers receive service before low priority ones. They proposed a mixed integer linear programming model and a metaheuristic approach based on adaptive large neighborhood search algorithm, used for small-scale and large-scale problems, respectively. Experimental results demonstrated excellent performance in both solution quality and computation time for the proposed algorithm. Li et al. (2022), taking into consideration, the acceleration of population aging and the limited availability of home care personnel, prioritized providing services to customers with more severe health conditions when planning routes. They developed a fuzzy multi-objective optimization model that maximizes the total priority of visiting customers while minimizing the total service cost. Additionally, they proposed a discrete multi-objective gray wolf optimizer, which effectively addressed the home care routing and scheduling problem. Zhang et al. (2021) established a dynamic VRP for prioritizing the delivery of “compensation for overtime” orders based on an Online To Offline (O2O) take-out platform and designed an improved iterative local search algorithm to solve the model. Wang and Huang (2022) identified priority distribution points through customer value assessment and ranking. They constructed a multi-objective cold chain logistics distribution path model and utilized a GA to effectively reduce distribution costs and enhance customer satisfaction. Li and Han (2023) used a cloud model to divide customers into three levels. They proposed an unmanned vehicle delivery paths optimization model based on constantly changing road conditions. The model was solved using a hybrid genetic-simulated annealing algorithm, with experimental results demonstrating its significant effectiveness.

In addition, VRP involves making other decisions (Mor and Speranza, 2022) and addressing challenges such as considering specific factors. For example, in order to reduce environmental impact, researchers attempted to reduce greenhouse gas emissions from vehicles and promote sustainable development by considering carbon emissions in route planning. They have established mathematical models that take carbon emissions into account and have studied optimization problems using algorithms such as PSO and ACO (Li et al., 2019; Cai et al., 2021; Chen et al., 2021; Guo et al., 2022).

Compare with the research on the pharmaceutical distribution routing problem and VRP, there is relatively little literature concentrates on the pharmaceutical distribution routing problem considering customer priority and carbon emissions, that is, the customer priority and carbon emissions minimization are simultaneously considered in the pharmaceutical distribution routing problem. The pharmaceutical distribution routing problem considering customer priority and carbon emissions is very important, the reasons are as follows: (1) Different types of customers have different levels of importance, considering customer priority will help to allocate the resources more reasonably and thus enhance the customer satisfaction. (2) Reducing carbon emissions is necessary for the logistics industry to achieve sustainable and green development.

GA is a population-based metaheuristic that has been successfully applied to solve a variety of combinatorial optimization problems, such as VRP and job shop scheduling problems (Chien and Lan, 2021; Wang et al., 2023; Liu et al., 2024). In order to further enhance intensification abilities of GA, some researchers have integrated local search into the GA, resulting in hybrid genetic algorithm (HGA). The HGA has been utilized to solve various NP-hard combinatorial optimization problems (Sun et al., 2023; Mahmoudinazlou and Kwon, 2024).

This article investigates the pharmaceutical distribution routing problem considering customer priority and carbon emissions. More specifically, the customer priority and carbon emissions minimization are simultaneously considered in the pharmaceutical distribution routing problem. The main contributions of this article are as follows:

1. (1) A mathematical model is established for the pharmaceutical distribution routing problem considering customer priority and carbon emissions. The objectives of the mathematical model are minimizing fixed cost, refrigeration cost, fuel consumption cost, carbon emission cost, and penalty cost for violating time windows.

2. (2) An effective HGA is proposed to solve the aforementioned problem, where an efficient local search based on VNS is developed and incorporated to enhance the intensification abilities.

The rest of the article is organized as follows: Section 2 states the pharmaceutical distribution routing problem considering customer priority and carbon emissions, and further establishes corresponding mathematical model. Section 3 designs a HGA to solve the abovementioned problem. Section 4 presents the experimental results, where the proposed HGA is compared with four optimization algorithms. Finally, Section 5 provides the conclusion and future work.

2. Mathematical model

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., 2021).

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., 2021). 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., 2021). 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. The proposed hybrid genetic algorithm (HGA)

This section introduces the proposed HGA for optimizing pharmaceutical distribution, taking into account customer priority and carbon emissions. The main components of the HGA are sequentially introduced, including encoding and decoding methods, population initialization, select operator, crossover based on adaptive probability, mutation based on adaptive probability, local search based on VNS, and recombination operator. Finally, the framework of the proposed HGA is presented.

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., 2023):

(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., 2023):

(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., 2010). Researchers have employed this method to successfully address numerous combinatorial optimization problems (Peng et al., 2019, 2023; Sun et al., 2023). 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. Experimental comparison results

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 (2021), 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, 2021), GA, improved ant colony algorithm (IACA; Guo et al., 2023), 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.

5. Conclusion and future work

This article pays great attention on factors such as customer priority and carbon emissions, and investigates the pharmaceutical distribution routing problem considering customer priority and carbon emissions. We establish a mathematical model for the problem, where the objectives are to minimize the carbon emissions and other optimization goals. This model aims to meet diverse customer demands and control enterprise operational costs, concurrently reducing carbon emissions incurred during pharmaceutical distribution processes. To solve this optimization problem, we propose an effective HGA. To enhance the algorithm performance of the proposed HGA, an efficient local search based on VNS is devised and incorporated into the HGA. Meanwhile, crossover with adaptive probability, and mutation with adaptive probability are used in the HGA. Finally, the effectiveness of the model and proposed HGA are validated through experiments.

In the future work, we will design more effective population initialization and local search strategies to further enhance the effectiveness of the proposed HGA. Moreover, from the perspective of multi-objective optimization, we will investigate the multi-objective pharmaceutical distribution routing problem considering customer priority and carbon emissions.