Multiobjective Route Optimization for Multimodal Cold Chain Networks Considering Carbon Emissions and Food Waste

Abstract

The cold chain logistics industry faces significant challenges in terms of transportation costs and carbon emissions. It is imperative to plan multimodal transportation routes efficiently to address these issues, minimize food waste, and reduce carbon emissions. This paper focuses on four key optimization objectives for multimodal cold chain transport: minimizing total transportation time, costs, carbon emissions, and food waste. To tackle these objectives, we propose a high-dimensional multiobjective route optimization model for multimodal cold chain networks. Our approach involves the development of a multiobjective evolutionary algorithm, utilizing Monte Carlo simulation and a one-by-one selection strategy. We evaluate the proposed algorithm’s performance by analyzing various convergence and distribution indicators. The average values for the minimum total transportation time, transportation cost, carbon emission cost, and cargo loss rate derived from the proposed algorithm ultimately converge to 6721.7, 5184.4, 301.5, and 0.21, respectively, demonstrating the effectiveness of the algorithmic solution. Additionally, we benchmark our algorithm against the existing literature to showcase its efficiency in solving high-dimensional multi-objective route optimization problems. Furthermore, we investigate the impact of different parameters, such as carbon tax rates, temperature, and cargo activation energy, on carbon emissions, and food waste. Moreover, we conduct a real-world case study to apply our approach to solving a practical business problem related to multimodal cold chain transportation. The insights gained from this research offer valuable decision-making support for multimodal carriers in developing low-carbon and environmentally friendly transportation strategies to efficiently transport perishable goods.

1. Introduction

As society and the global economy evolve, there has been a marked increase in health awareness, leading to a stronger emphasis on the freshness and quality of food. In recent years, the global production of fresh agricultural products has steadily risen, with projections suggesting a potential 60% increase by 2050. This growth is driven by improving living standards and increasing awareness of health among consumers, which has driven the expansion of supply chains and the proliferation of global retail networks [1]. Technological advances, including the Internet of Things (IoT), artificial intelligence (AI), and digital twins, have also played a pivotal role in improving the monitoring and tracking of product transport [2]. Moreover, the rise of fresh e-commerce platforms, vegetable and fruit delivery services, and online shopping has significantly boosted the fresh e-commerce market, further fueling the demand for efficient cold chain logistics. As a result, the trade of perishable goods has emerged as a critical component of international commerce.

According to reports from the Food and Agriculture Organization, more than 14% of food is lost during transport from production sites to retailers, and approximately 30% of fresh food spoils in transit due to inadequate temperature regulation [3,4]. The perishable nature of fresh goods requires the use of specialized personnel, equipment, and technologies to continuously monitor environmental conditions, cargo quality, and storage throughout the transportation process. Ensuring the integrity, safety, and freshness of these goods not only enhances customer satisfaction but also strengthens the competitiveness of logistics providers [5]. However, compared to conventional goods, transporting fresh products is more time consuming and costly. The need for temperature-controlled environments to maintain product quality further increases carbon emissions, which conflicts with lowcarbon environmental goals and the imperatives of a sustainable economy.

The transportation sector is a major contributor to energy consumption and carbon dioxide emissions, posing significant challenges to achieving carbon neutrality [6]. In pursuit of a low-carbon transformation in transportation, upgrading infrastructure and promoting multimodal transportation have emerged as sustainable solutions. These approaches reduce transportation costs and times, enhance logistics standards, alleviate road congestion, and minimize environmental impacts. Driven by supportive government policies and market demands, the cold chain logistics industry is poised to adopt multimodal transportation practices. This is a composite transportation model, combining multiple modes of transport that connect and transfer seamlessly to complete the entire cargo transportation process. By optimizing resource allocation, improving organizational structures, and reducing dependence on traditional road and air transport, multimodal cold chain logistics can offer more efficient distribution routes. This approach not only mitigates product loss and reduces carbon emissions but also enhances overall logistics efficiency.

Existing research underscores the critical importance of managing food waste in cold chain logistics. Shashi et al. [7], in their systematic review of 1189 studies on the Food Cold Chain (FCC), identified food waste management as one of the most frequently cited topics in FCC literature. Wu and Hsiao [8], using a failure mode and effects analysis (FMEA), identified improper temperature control as the leading cause of food quality degradation in cold chain systems. Multimodal temperature-controlled transportation, highlighted by Behdani et al. [9], offers an economical and efficient alternative to single-mode road transport, addressing key challenges in managing perishable goods. Zhang et al. [10] proposed a value-based management (VBM) approach within a decision-making framework for multimodal cold chains, optimizing transport mode selection and scheduling by leveraging the characteristics of perishable goods and market conditions.

The inherent characteristics of perishable goods, such as order volume, temperature control requirements, and waste cost coefficients, introduce complexity into distribution scheduling and delivery route optimization models. Zhang et al. [11] addressed a vehicle routing problem (VRP) for multiproduct frozen food distribution by utilizing a genetic algorithm (GA) that accounts for unit volume and perishability coefficients. Chen et al. [12] investigated the scheduling of distribution and the optimization of route for cold food chains under increasing demand, taking into account both the volume of orders and the temperature control requirements. Zheng et al. [13] integrated food loss costs and freshness satisfaction into a multimodal transport mode selection model to optimize cold chain food logistics for efficient delivery. Urban logistics systems add another layer of complexity to cold chain VRP models. Chaofan et al. [14] developed a VRP for urban cold chain distribution that incorporates simultaneous delivery and pickup to minimize total transportation costs using a GA to solve the model. Franceschetti et al. [15] examined pollution-aware VRP models, proposing a metaheuristic algorithm based on adaptive large neighborhood search (ALNS) to address time-dependent pollution in urban traffic. Ma et al. [16] combined order selection and a time-dependent VRP with time windows for perishable goods delivery, solving the model with a hybrid ant colony optimization algorithm. Liu [17] extended this work with a time-dependent cold chain multimodal transportation route optimization model, which addresses both static and dynamic realtime optimization. Transportation planning models for cold chain products are typically shaped by preservation constraints or penalty costs, focusing on a single aspect. However, most literature lacks consideration of the temperature variations in quality control modeling, leading to insufficiently comprehensive problem formulations.

Incorporating multiple objectives in cold chain optimization is essential to ensure efficiency and sustainability, as single-objective models often fall short. Existing studies commonly convert multiobjective optimization into a single-objective problem by applying weights to each objective or using chance constraints. However, these methods can lead to a loss of feasible solutions and restrict the decision makers’ range of choices. Roghanian and Cheraghalipour [18] tackled this challenge by developing a closed-loop citrus supply chain optimization model with objectives that minimize cost, maximize demand responsiveness, and reduce carbon emissions. They employed five metaheuristic algorithms, with the multiobjective tree growth algorithm (MOTGA) outperforming others. Tirkolaee et al. [19] also recognized the importance of multiobjective optimization, introducing a pollution routing problem with cross-docking using a biobjective mixed-integer linear programming model, which as solved with MOSA and NSGA-II algorithms. Zulvia et al. [20] proposed a multi-objective green VRP for perishable products, incorporating operational costs, deterioration costs, carbon emissions, and customer service levels, and they developed a many-objective gradient evolution algorithm (MOGE), which outperformed other algorithms in terms of diversity and convergence. Wu et al. [21] considered multiple factors in their model, including transportation costs, penalties, overloading costs, carbon tax costs, and customer satisfaction. They formulated the objective function as the ratio of total costs to customer satisfaction and applied an enhanced A* algorithm alongside an ant colony algorithm to solve the model. The results verified the model’s effectiveness, efficiency, and accuracy. Zhang et al. [22] developed a low-carbon cold chain logistics route optimization model tailored to China’s traffic conditions, factoring in the combined impacts of traffic conditions on total distribution costs, product freshness, and carbon emissions. To address the complexities of this multidimensional discrete model, they proposed an improved discrete firefly algorithm (IDFA). Similarly, Liao et al. [23] designed an urban cold chain logistics network with a control strategy aimed at balancing total cost and customer satisfaction, thereby reducing the risk of outbreaks associated with cold chain pollutants. To solve this model, they introduced a specially initialized multiobjective evolutionary algorithm (MOEA-SI). It is evident that the research on optimization algorithms for discrete multiobjective problems, typically involving two or three objectives, has been widely explored in the literature. However, there is comparatively less focus on optimization problems with four or more objectives in this field.

Our study addresses the research gaps by developing a high-dimensional multiobjective route optimization model for multimodal cold chain transportation networks. The model aims to minimize total transportation costs, time, carbon emissions, and food loss rates. Transportation costs include on-route, trans-shipment, and storage costs. To solve the model, we propose a multiobjective evolutionary algorithm based on Monte Carlo simulation and a one-by-one selection strategy (MC-ObOEA), which is designed to handle the uncertainty in transportation time. The effectiveness of the algorithm will be demonstrated through a case study using real-world business data. The main contributions of this paper are as follows:

• The model accounts for the accumulation of total transportation time and temperature variations during transit to highlight the differences between cold chain cargo transportation and traditional cargo transportation. The Weibull function is used to model the changing cargo loss rate due to temperature fluctuations.
• To make the model more realistic, random variables are incorporated to represent transportation time uncertainties, and Monte Carlo simulation is employed for the effective treatment of these uncertainties. Additionally, the shift periods of transportation modes are explicitly modeled.
• A carbon pricing function is introduced to convert carbon emissions into carbon emission costs, thereby integrating environmental impact into the transportation modeling process.
• To solve the high-dimensional, multiobjective optimization problem, we propose the MC-ObOEA algorithm. The algorithm evaluates and selects individuals based on convergence and distribution metrics, ensuring that the Pareto front solution set with superior performance is obtained.

The remainder of the paper is organized as follows: Section 2 provides a formal description of the problem. Section 3 introduces the high-dimensional MC-ObOEA algorithm. Computational experiments are discussed in Section 4, followed by a case study in Section 5. Finally, Section 6 presents the conclusions of the study.

2. A Multiobjective Multimodal Cold Chain Transport Route Optimization Model

2.1. Problem Description

The multimodal transportation network analyzed in our study encompasses highway, railway, and waterway transport modes. The decision maker may select one or multiple transportation modes to deliver goods from the distribution center to the final destination. Each mode significantly impacts transportation time, costs, carbon emissions, and food losses in cold chain logistics. Highway transport is subject to variability caused by road congestion and traffic accidents, while waterway transport is influenced by factors such as extreme weather conditions, port congestion, and ship collisions, all of which introduce uncertainties in transportation time. Both the waterways and the railway modes follow fixed daily schedules, which affect not only the timing of subsequent departures but also limit cargo waiting times at transshipment nodes, leading to increased holding costs.

In the context of carbon reduction policies and efforts to promote a low-carbon transformation in the transportation sector, it is essential to account for the cost of carbon emissions during the transportation process. Moreover, in the supply chain of fresh and perishable goods, switching between transport modes often extends transit times and accelerates perishable product degradation, resulting in significant food waste. Taking these factors into account, this study performs a single-source, single-destination route optimization with the objectives of minimizing total transportation time, overall transportation costs, carbon emissions, and food loss rates. The study is based on the following assumptions:

1.

Transported cargoes are indivisible and must remain intact during delivery.

2.

The time and cost of transitioning between different transport modes are known.

3.

The arrival time of the cargoes at the nodes marks the beginning of loading, unloading, and trans-shipment.

4.

After loading, unloading, and trans-shipment, the subsequent transport leg begins according to the nearest available schedule of the chosen mode.

5.

The transport time distributions for each stage of the journey are known.

6.

Waiting times during transportation are solely due to constraints of the transport mode shift.

7.

The capacities, facilities, equipment, and personnel of the transport nodes meet the requirements for feasible transshipment.

8.

The temperature of the goods remains constant while being transported by a single mode of transport; however, during transit, the temperature may change to some extent.

2.2. Parameter Definition

We represent the multimodal stochastic transportation network as a graph $G=(N,\tilde{E},M)$, where N is the set of nodes, with nodes o and d representing the nodes of origin and destination, respectively. M denotes the set of available transport modes, and $\tilde{E}$ is the set of arcs, which is defined as $\tilde{E}=\left\{{\tilde{e}}_{i,j}^a|i,j\in N,a\in M\right\}$. Each arc ${\tilde{e}}_{i,j}^a$ is characterized by a set of parameters ${\tilde{t}}_{i,j}^a$, $c_{i,j}^a$, $c_{e_{i,j}}^a$, and $D_{i,j}^a$, where the transport times ${\tilde{t}}_{i,j}^a$ follow a specific probability distribution. Detailed definitions of the parameters are provided below:

$q:$

Cargo volume.

$p:$

Carbon tax coefficient.

$c_{i,j}^a:$

Transportation cost from node i to node j with transport mode a.

${\tilde{t}}_{i,j}^a:$

Transportation time from node i to node j with transport mode a.

${\delta }_{i,j}^a:$

Carbon emissions coefficient from node i to node j with transport mode a.

$v_{i,j}^a:$

Average speed from node i to node j by transport mode a.

$T_i^{dep}:$

Cargo departure time at node i (assuming that the departure time at the origin $T_o^G$ is given).

$T_i^{arr}:$

Cargo arrival time at node i.

$T_i^{fin}:$

Cargo trans-shipment completion time at node i.

$c_i^{a,b}:$

Cost of changing the transport mode from a to b at node i.

${\tilde{t}}_i^{a,b}:$

Time to change the transport mode from a to b at node i.

${\xi }_i^{a,b}:$

Carbon emissions coefficient of change in transport mode from a to b at node i.

${\Gamma }_{i,j}^a:$

Gives the schedule for using the transport mode a from node i to node j, where ${\Gamma }_{i,j}^a=\left(\dots ,{\mathrm{\Phi }}_{n-1},{\mathrm{\Phi }}_n,\dots \right)$ and ${\mathrm{\Phi }}_n$ denote the departure time of the n-th change.

$t_i^{del}:$

Cargo transshipment time at node i.

$f_i:$

The correlation between cargo holding costs and trans-shipment time is expressed as a piecewise function, where the cost varies depending on the duration of the holding period. (For example, the holding cost is set at 0.01 CNY per minute for the first three hours, increasing to 0.02 CNY per minute for three to five hours, with different rates applied beyond this time.) These rates also differ depending on the transport mode used.

A set of binary decision variables is defined as follows:

$$\begin{array}{cc}\hfill & x_{i,j}^a=\left\{\begin{array}{c}1,\mathrm{if}\mathrm{the}\mathrm{cargo}\mathrm{is}\mathrm{shipped}\mathrm{from}\mathrm{node}i\mathrm{to}\mathrm{node}j\mathrm{via}\mathrm{the}\mathrm{transport}\mathrm{mode}a\hfill \\ 0,\mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & y_j^{a,b}=\left\{\begin{array}{c}1,\mathrm{if}\mathrm{the}\mathrm{cargo}\mathrm{is}\mathrm{transferred}\mathrm{from}\mathrm{transport}\mathrm{mode}a\mathrm{to}\mathrm{mode}b\mathrm{at}\mathrm{node}j\hfill \\ 0,\mathrm{otherwise}\hfill \end{array}\right.\hfill \end{array}$$

(2)

2.3. Optimization Model Formulation

When the cargo reaches node i for trans-shipment, it departs at the latest scheduled departure time of the selected transport mode after the completion of the transfer process. The departure time of the cargo is determined as follows:

$${\mathrm{\Phi }}_n=\phi \left(T_i^{fin},{\Gamma }_{i,j}^a\right),{\mathrm{\Phi }}_{n-1}<T_i^{fin}\le {\mathrm{\Phi }}_n$$

(3)

The moment when the transition is completed after the shipment arrives at a node i is as follows:

$$T_i^{fin}=(T_i^{arr}+{\tilde{t}}_i^{a,b})y_i^{a,b},\forall \left(i\right)\in N\setminus \left\{o\right\}$$

(4)

Then, the departure moment $T_i^{dep}$ at node i is as follows:

$$T_i^{dep}=\sum _{b\in M}\sum _{j\in N\setminus \left\{o\right\}}\phi \left(T_i^{fin},{\Gamma }_{i,j}^{\mathrm{b}}\right)x_{i,j}^a,\forall \left(i\right)\in N\setminus \left\{d\right\}$$

(5)

The moment $T_j^{arr}$ at which the shipment arrives at node j is shown below:

$$T_j^{arr}=(T_i^{dep}+{\tilde{t}}_{i,j}^a)x_{i,j}^a,\forall \left(i\right)\in N\setminus \left\{o\right\}$$

(6)

The detention time $t_i^{del}$ of the cargo at node i is as follows:

$$t_i^{del}=T_i^{dep}-T_i^{fin},\forall \left(i\right)\in N\setminus \{o,d\}$$

(7)

Therefore, the total transportation time T can be calculated in the following:

$$T=\sum _{i\in N}\sum _{j\in N\setminus \left\{i\right\}}\sum _{a\in M}{\tilde{t}}_{i,j}^a{x}_{i,j}^a+\sum _{i\in N\setminus \{o,d\}}\sum _{a\in M}\sum _{b\in M}{\tilde{t}}_i^{a,b}{y}_i^{a,b}+\sum _{i\in N\setminus \{o,d\}}\sum _{a\in M}\sum _{b\in M}{t}_i^{del}{y}_i^{a,b}$$

(8)

Transportation costs are calculated based on three primary components: in-transit road transportation costs, transshipment costs, and additional storage costs incurred from holding cargo at trans-shipment nodes due to scheduling constraints.

$$C=\sum _{i\in N}\sum _{j\in N\setminus \left\{i\right\}}\sum _{a\in M}{c}_{i,j}^a{x}_{i,j}^a+\sum _{i\in N\setminus \{o,d\}}\sum _{a\in M}\sum _{b\in M}{c}_i^{a,b}{y}_i^{a,b}+\sum _{i\in N\setminus \{o,d\}}\sum _{a\in M}\sum _{b\in M}{f}_i{t}_i^{del}{y}_i^{a,b}$$

(9)

During the model construction, a carbon tax function is incorporated to translate carbon emissions into associated costs. Carbon emissions are primarily derived from two sources: transportation and trans-shipment. Taking both into account, the total carbon emission cost, $C_e$, is calculated using the following formula:

$$C_e=(\sum _{i\in N}\sum _{j\in N\setminus \left\{i\right\}}\sum _{a\in M}{\tilde{t}}_{i,j}^a{v}_{i,j}^a{\delta }_{i,j}^a{x}_{i,j}^a+\sum _{i\in N\setminus \{o,d\}}\sum _{a\in M}\sum _{b\in M}{\tilde{t}}_i^{a,b}{\xi }_i^{a,b}{y}_i^{a,b})qp$$

(10)

The key distinction between cold chain transport and traditional cargo transport lies in the gradual decline of product quality during storage and transit. The percentage of damaged goods reflects the loss of cargo due to quality degradation caused by the cumulative effects of transportation time and temperature fluctuations during transit and transshipment. The rate of deterioration due to temperature changes can be modeled using the Arrhenius equation: $K=A·exp(-E\alpha /RF)$, where K is the reaction rate constant, A is the Arrhenius constant (or frequency factor), $E\alpha$ represents the activation energy, and R is the gas constant. The reaction rates and temperatures of the cargo during transport and transit are denoted by $k_1$, $F_1$, and $k_2$, $F_2$, respectively. $k_1$ and $k_2$ are calculated using Equation (11).

$$k_1=A·exp\left[{\displaystyle \frac{-E_{\alpha }}{R\left(F_1+273.15\right)}}\right],k_2=A·exp\left[{\displaystyle \frac{-E_{\alpha }}{R\left(F_2+273.15\right)}}\right].$$

(11)

The relationship between $F_1$, $F_2$, and the temperature variation during transit $\Delta F$ is expressed in Equation (12).

$$\Delta F=\left|F_1-F_2\right|$$

(12)

Our cargo loss rate model is based on the time-dependent relationship of the cargo loss rate, which is represented by a Weibull function introduced by Giri et al. [24]. This relationship is expressed in Equation (13):

$$D_{\left(t\right)}=\left\{\begin{array}{c}1-e^{-kt},\mathrm{if}t>0\hfill \\ 0,\mathrm{else}\hfill \end{array}\right.$$

(13)

Therefore, the cumulative cargo loss rate during transportation can be calculated using the following formula:

$$D=\sum _{i=1}^n\sum _{j=1}^n\sum _{a=1}^M\left[1-exp(-k_1{\tilde{t}}_{i,j}^a)x_{i,j}^a+(1-exp(-k_2({\tilde{t}}_i^{a,b}+t_i^{del})))exp(-k_1{\tilde{t}}_{i,j}^a)y_i^{a,b}\right]$$

(14)

The optimization model is designed with four primary objectives: minimizing total transportation time, reducing total transportation costs, minimizing carbon emission costs, and minimizing cargo loss rates. These objectives are incorporated into a high-dimensional, multiobjective, multimodal transport route optimization model, while the previously discussed constraints are implicitly included but not restated in the following formulation.

Objectives:

$$minT,minC,min{C}_e,\mathrm{and}minD$$

(15)

Constraints:

$$\begin{array}{cc}\hfill & \sum _{j\in N\setminus \left\{i\right\}}\sum _{a\in M}{x}_{i,j}^a\le 1,\forall i\in N\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & \sum _{i\in N\setminus \left\{o\right\}}\sum _{a\in M}{x}_{o,i}^a=\sum _{i\in N\setminus \left\{d\right\}}\sum _{a\in M}{x}_{i,d}^a=1\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill & \sum _{i\in N\setminus \left\{o\right\}}\sum _{a\in M}{x}_{i,o}^a=\sum _{i\in N\setminus \left\{d\right\}}\sum _{a\in M}{x}_{d,i}^a=0\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill & \sum _{a\in M}\sum _{b\in M}{y}_i^{a,b}\le 1,\forall \left(i\right)\in N\setminus \{o,d\}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill & \sum _{i\in N\setminus \left\{j\right\}}\sum _{a\in M}{x}_{i,j}^a=\sum _{i\in N\setminus \left\{j\right\}}\sum _{a\in M}{x}_{j,i}^a,\forall j\in N\setminus \left\{o,d\right\}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & \sum _{k\in N\setminus \left\{d\right\}}{x}_{k,i}^a+\sum _{j\in N\setminus \left\{o\right\}}{x}_{k,i}^b\ge 2y_i^{a,b},\forall \left(i\right)\in N\setminus \{o,d\},\forall a,b\in M\hfill \end{array}$$

(21)

Equation (16) guarantees that only a single-cargo departure occurs at each node. Equations (17) and (18) ensure that cargo departs from the origin node o (the initial shipment point) and eventually reaches the destination node d (the final delivery point). Equations (19) and (20) stipulate that the cargo is loaded and unloaded at one transit node at a time. Lastly, Equation (21) ensures that the flow balance is maintained across all nodes.

3. A High-Dimensional Multiobjective Multimodal Route Optimization Algorithm

We introduce a high-dimensional, multiobjective, multimodal transportation route optimization algorithm, called the Multiobjective Evolutionary Algorithm Based on Monte Carlo Simulation and One-by-One Selection Strategy (MC-ObOEA), to solve our proposed model. The algorithm is built on the core framework of a genetic algorithm [25] and incorporates the Monte Carlo (MC) simulation [26] to account for uncertainties in transportation time. The algorithm evaluates individuals using both convergence and distribution indicators. During the offspring generation process, individuals are selected primarily based on convergence indicators, and for each selected individual, similar individuals are eliminated using distribution indicators to maintain a balance between convergence and diversity in the population. To further enhance the performance of the solution set, a boundary preservation mechanism based on angle solutions is also incorporated, improving the extension and coverage of the results.

3.1. General Algorithm Framework

Algorithm 1 outlines the overall framework, which consists of the following steps. First, a randomly initialized population is generated based on the initial transport network and arc segment data. The offsprings are then produced using cross and mutation operators. The parent and offspring populations are then merged to form a combined population, and the fitness of each individual is evaluated. Following this, the convergence and distribution metrics are computed across the population. Finally, a one-by-one selection process is applied to form the new population. Two critical components of this framework are the calculation of convergence and distribution metrics (lines 7–9), as well as the selection strategy (line 10), both of which are detailed in the following sections.

Algorithm 1 MC-ObOEA framework.

Input: P (initial population), N (population size)  1: $t=0$  2: Initializing the population $P\left(t\right)$;  3: while Failure to meet termination conditions do  4:    Generation of offspring populations $P^{\prime }\left(t\right)$ using crossover, mutation operators;  5:    $Q\left(t\right)=P\left(t\right)\cup P^{\prime }\left(t\right)$;  6:    Calculate the multi-objective adaptation value for individuals in $Q\left(t\right)$;  7:    for all $x_i\in Q\left(t\right),i=1,\dots ,2N$ do  8:      Compute the convergence metric $c\left(x_i\right)$ and the distribution metric $d\left(x_i\right)$;  9:   end for 10:   $P(t+1)$ selects $\left(Q\right(t),c,d)$ one by one; 11:   $t=t+1$ 12: end while 13: Choose the unoccupied solution in $P\left(t\right)$ to form $FS$ 14: Output $FS$

3.2. Convergence and Distribution Metrics

3.2.1. Convergence Metric

In calculating the convergence metric, all optimization objectives are treated with equal importance and are combined into a single scalar value. The general formulation of the convergence metric is presented as follows:

$$c\left(x\right)=agg[f_1\left(x\right),\dots ,f_M\left(x\right)]$$

(22)

Three representative methods are used to calculate the convergence metric. The first method involves determining the Chebyshev distance from each individual to the globally optimal solution, which is denoted as $CdI$.

$$c\left(x\right)=\underset{1\le m\le M}{max}{f}_m\left(x\right)-z_m^{\ast }$$

(23)

where $z^{\ast }={\left\{z_1^{\ast },\dots ,z_M^{\ast }\right\}}^T$ is the globally optimal value, and $z_m^{\ast }=\underset{x_i\in Q\left(T\right)}{min}{f}_m\left(x_i\right)$. The second method calculates the Euclidean distance between each individual and the global optimal value, which is denoted as $EdI$.

$$c\left(x\right)=\sqrt{{\sum }_{m=1}^M{\left[f_m\left(x\right)-z_m^{\ast }\right]}^2}$$

(24)

The third method calculates the Euclidean distance between each individual and the global base value, which is referred to as $EdN$.

$$c\left(x\right)=1/\phantom{1\sqrt{{\sum }_{m=1}^M{\left[f_m\left(x\right)-z_m^{nad}\right]}^2}}\sqrt{{\sum }_{m=1}^M{\left[f_m\left(x\right)-z_m^{nad}\right]}^2}$$

(25)

where $z^{nad}={\left\{z_1^{nad},\dots ,z_M^{nad}\right\}}^T$ is the global base value, and $z_m^{nad}=\underset{x_i\in Q\left(t\right)}{max}{f}_m\left(x_i\right)$. Since our model is a minimization problem, the inverse of the distance is used as an indicator of convergence.

Figure 1 illustrates the contour plots formed by different convergence metrics in a two-dimensional objective space $(f_1,f_2)$. Individuals located on the same contour line share identical convergence metric values. When the contours of the chosen convergence metric align closely with the shape of the Pareto frontier of the optimization problem, the resulting solution set tends to be of higher quality.

3.2.2. Distribution Metric

The distribution metric is represented as a vector where each element denotes the distance between one individual and every other individual. The cosine similarity, calculated as the cosine of the angle between two vectors, is used to assess their similarity. A cosine similarity of 1 indicates identical vector directions, while a value of 0 signifies perpendicular vectors. The range of cosine similarity spans [0, 1], making it particularly suitable for comparing individuals in high-dimensional spaces. Figure 2 illustrates a step by-step selection process using cosine similarity as a distribution metric, with $cos\left(\theta \right)$ as the distribution threshold. In the example, the individual A is selected, leading to the elimination of C. Then, the individuals B and D are selected, and E is cleared. Finally, F is selected, resulting in the removal of I, H, and G. This method effectively reduces the number of solutions that exhibit dominant resistance, such as I, H, and G.

Therefore, the following distributive indicators are introduced.

$$d\left(x_i\right)=\left\{d_1\left(x\right),\dots ,d_{\left|Q\left(t\right)\right|}\left(x_i\right)\right\},i=1,\dots ,\left|Q\left(t\right)\right|$$

(26)

where

$$d_j\left(x_i\right)=1-cos\left({\theta }_{ij}\right),j=1,\dots ,\left|Q\left(t\right)\right|$$

(27)

represents the distance between $x_i$ and $x_j$, and

$$cos\left({\theta }_{ij}\right)={\displaystyle \frac{{\sum }_{m=1}^M\left[f_m\left(x_i\right)-z_m^{\ast }\right]\times \left[f_m\left(x_j\right)-z_m^{\ast }\right]}{\sqrt{{\sum }_{m=1}^M{\left[\left(x_i-z_m^{\ast }\right)\right]}^2\times \sqrt{{\sum }_{m=1}^M{\left[f_m\left(x_i\right)-z_m^{\ast }\right]}^2}}}}.$$

(28)

It is straightforward to observe that $d_j\left(x_i\right)$ takes values within the range [0, 1], where smaller values of $d_j\left(x_i\right)$ indicate that $x_i$ is closer to $x_j$. Furthermore, $d_j\left(x_i\right)$ is equal to $d_i\left(x_j\right)$, and $d_i\left(x_i\right)$ is always 0. During the selection process, once an individual is selected, all individuals in the population with $d_j\left(x_i\right)$ below a predefined threshold are removed. If a distribution metric based on Euclidean distance is employed, the formula for calculating the distance between $x_i$ and $x_j$ is given as follows:

$$d_j\left(x_i\right)=\sqrt{{\sum _{m=1}^M\left[f_m\left(x_i\right)-f_m\left(x_j\right)\right]}^2},j=1,\dots ,\left|Q\left(t\right)\right|.$$

(29)

3.3. Selection Strategy

The flowchart of the MC-ObOEA algorithm is illustrated in Figure 3. The gray boxes detail the process of individual selection, carefully outlining the steps of the selection strategy. Using the convergence and distribution metrics introduced in Section 3.2, the individuals in population $Q\left(t\right)$ are categorized into three distinct sets: (1) the preselected solution set $Q_s$, (2) the set of solutions $Q_{th}$ eliminated by distribution thresholds, and (3) the set of dominated solutions $Q_d$. Initially, a boundary preservation mechanism selects corner individuals, denoted as $x_m^{coner}$, for inclusion in the $Q_s$ solution set. In an optimization problem with M objective functions, a corner solution is defined as the one that simultaneously minimizes $k_c$ objectives, where $k_c<M$. We define a corner solution (or marginal individual) as an individual with the smallest value of the scalar function in the current population of the objective functions $k_c$. These corner individuals are identified using M scalar functions, which replace convergence metrics in the selection process.

$$c_m\left(x\right)=agg\left[f_1\left(x\right),\dots ,f_{m+1}\left(x\right),\dots ,f_M\left(x\right)\right],m=1,\dots ,M$$

(30)

For example, when using $EdI$ as the convergence metric, the scalar function is given by

$$c_m\left(x\right)=\sqrt{{\sum }_{i=1,i\ne m}^M{\left[f_i\left(x\right)-z_i^{\ast }\right]}^2},m=1,\dots ,M.$$

(31)

Therefore, the marginal individual $x_m^{coner}$ can be obtained as follows:

$$x_m^{coner}=arg\underset{x_i\in Q\left(t\right)}{min}{c}_m\left(x_i\right),m=1,\dots ,M$$

(32)

Prioritizing the selection of these edge individuals significantly improves the algorithm’s ability to explore the entire Pareto frontier, as these individuals are inherently non-dominated. Afterward, the remaining individuals of the population $Q_t$ are selected for inclusion in $Q_s$ one by one. In particular, when an individual $x_{unique}$ is added to the $Q_s$ solution set, any individual that is dominated by or too similar to it is removed from the population using the removal operator $\zeta$. This ensures that dominated individuals are not selected over non-dominated ones. The selection and removal operations continue until $Q_t$ is fully depleted. If, after this initial depletion, the number of individuals in $Q_s$ exceeds the population size N, only the first N individuals are retained for population $P\left(t\right)$. Otherwise, the remaining individuals in $Q_{th}$ and $Q_d$ will compete for selection in $Q_s$ again. During this process, once $Q_t$ is empty, the individuals of $Q_{th}$ and $Q_d$ are reintroduced into $Q_t$ for another selection round. Finally, based on the number of individuals in $Q_s$, the distribution threshold is updated as follows:

$${\zeta }^{t+1}={\zeta }^t\times e^{{\displaystyle \frac{ratio/\phantom{ratioRATIO-1}RATIO-1}{M}}},$$

(33)

where $ratio=\left|Q_s\right|/\phantom{\left|Q_s\right|N}N$ is the ratio of the number of preselected solutions to the population size. In our study, ${\zeta }^1$ is initialized to 1, and the threshold $RATIO\in \left[0,2\right]$ is used to control the number of preselected solutions. When $RATIO=1$, it provides the optimal balance between convergence and distribution performance. Using the adaptive adjustment method in Equation (33), if the number of preselected solutions falls below $R\times N$, ${\zeta }^t$ is reduced; conversely, if it exceeds this value, ${\zeta }^t$ is increased. This ensures that the number of preselected solutions remains approximately at $R\times N$. However, when the number of dominated individuals in $Q_t$ exceeds $(2-R)\times N$, it becomes impossible to select N individuals for $Q_s$, regardless of how ${\zeta }^t$ is adjusted. In such cases, the value of ${\zeta }^{t+1}$ remains unchanged.

4. Computational Experiments

4.1. Experimental Design and Parameters

This paper constructs a multimodal transportation network consisting of 15 nodes, where H represents highways, R represents railroads, and W represents waterways, as illustrated in Figure 4. Assuming that the transportation time for each segment follows a normal distribution, the specific data are provided in

Table 1. Mutimodal transportation network dataset.

oA;[25,3]891;[42,4]700;∼;	EF;[47,5]1650;∼;∼;
oB;[25,3]924;[42,3]700;∼;	EG;[47,5]891;∼;∼;
oC;[28,5]957;[70,6]372.75;∼;	EH;[31,3]1551;[43,3]980;[102,5]530.25;
AD;[48,5]1650;∼;∼;	EJ;∼;∼;[95,6]509.25;
AE;[25,4]924;[42,3]700;∼;	EM;∼;∼;[95,6]1120;
AI;∼;[38,5]1330;[100,6]556.5;	FI;∼;[38,5]1330;[100,6]551.25;
BA;[25,3]891;[43,4]700;∼;	FJ;[25,3]891;[43,3]700;∼;
BC;∼;[30,3]1120;∼;	FL;[25,3]891;[43,4]700;∼;
BD;[28,5]957;∼;[69,6]372.75;	GI;[25,3]924;[44,3]700;∼;
BE;[30,3]1551;[43,3]1015;[103,5]530.25;	GJ;∼;[30,3]1120;∼;
BF;[25,3]924;[42,4]700;∼;	GL;[24,3]924;[42,4]700;∼;
BG;[25,3]891;[42,4]700;∼;	HI;[31,3]1518;[43,3]1015;[102,5]530.25;
BH;[31,3]1551;[43,3]980;[103,5]540.75;	HJ;[28,5]957;∼;[70,6]372.75;
CD;[25,3]924;[43,4]700;∼;	HL;[24,3]924;[42,4]700;∼;
CE;[25,3]891;[43,4]700;∼;	IK;[25,4]891;[42,4]700;∼;
CJ;∼;[37,5]1295;[100,6]556.5;	IL;[30,3]1551;[43,3]101;[102,5]530.25;
DF;[25,3]924;[42,4]700;∼;	IM;∼;[37,4]1330;[101,5]551.25;
DG;[25,3]924;[43,3]700;∼;	JK;[25,3]11,891;[42,4]700;∼;
DH;[25,3]891;[42,4]700;∼;	JL;[30,4]1518;[43,3]1015;[103,4]530.25;
DI;∼;∼;[95,6]504;	JM;[47,4]1650;∼;∼;
DK;[25,3]891;[42,4]700;∼;	Kd;[30,4]1518;[43,3]1015;[102,5]535.5;
ML;[26,3]924;[43,4]700;∼;	LK;[47,4]1650;∼;∼;
Md;∼;[37,5]1295;[100,6]551.25;	Ld;[38,5]1330;[100,6]551.25;∼;

1 Data order: Origin and destination; road transport time and costs; rail transport time and costs; waterway transport time and costs. 2 Data format for randomly distributed parameters of transportation time: [average, standard deviation]. 3 Empty data: Represented as ∼.

. Additional data for nodal transit transitions, among other factors, are detailed in

Table 2. Carbon emissions per unit and average speed of various transport modes.

Transport Mode	Carbon Emissions (kgCO2/t km)	Average Speed (km/h)	Transport Mode
highway	0.044	80	highway
railway	0.0091	40	railway
waterway	0.0127	25	waterway

,

Table 3. Unit transition time between transport modes.

Transport Mode	Highway Transition Time (h)	Railroad Transition Time (h)	Waterway Transition Time (h)
highway	1	1.5	2
railway	1.5	2	4
waterway	2	4	3

,

Table 4. Unit transition cost between transport modes.

Transport Mode	Highway Transition Cost (CNY)	Railroad Transition Cost (CNY)	Waterway Transition Cost (CNY)
highway	30	30	40
railway	30	40	50
waterway	40	50	50

and

Table 5. Unit transition carbon emissions coefficients between transport modes.

Transport Mode	Highway Transition Coefficient kgCO2/t	Railroad Transition Coefficient kgCO2/t	Waterway Transition Coefficient kgCO2/t
highway	0.125	0.128	0.117
railway	0.128	0.115	0.113
waterway	0.117	0.113	0.12

. The interval between shifts is assumed to be consistent across all transport modes at each node. The MC-ObOEA algorithm is based on the genetic algorithm framework. The time complexity of the genetic algorithm can be approximated as $\mathrm{O}\left(N_d\times N\times F\right)$, where $N_d$ represents the number of decision variables, N is the population size, and F accounts for the time required to evaluate the fitness function and perform other operations with fixed time costs. The space complexity is approximately $\mathrm{O}\left(P\times S\right)$, where S is the space occupied by each individual. Given the more complex operations in MC-ObOEA, its overall complexity is higher, making the choice of algorithm parameters crucial. The algorithm and model parameters used in this study are based on the relevant literature [26,27], with the specific values presented in

Table 6. Parameters of the algorithm and model.

Parameter Definition	Initial Value	Parameter Definition	Initial Value
Network node size	15	Carbon tax coefficient p	30
Number of individuals in population	80	Temperature during transportation $F_1$ (°C)	5 °C
Number of iterations $Nd$	100	Variation in transit temperature range $\Delta F$(°C)	1 °C
Crossover probability $Pc$	0.8	Cargo activation energy $E\alpha$	34 KJ/mol
Mutation probability $Pm$	0.2	Cargo weight q	25 t
Simulation runs $MC$	500	Maximum reaction rate $K_{max}$	5 × 1014
$RATIO$	1	Gas constant R	8.314
Correlation coefficients between cargo storage costs and detention times			1

.

In addition to the 15-node network, 30-node, and 50-node networks were randomly generated following the same rules. The departure time for cargo from the origin is set to 7:30 a.m. In the MC-ObOEA algorithm—to achieve a balance between convergence and distribution performance—and the convergence indicator $EdI$ is utilized for optimization, while cosine similarity ($Cs$) is employed as the distribution indicator.

4.2. Impact of Convergence Indicators on the Performance of MC-ObOEA

In Equation (34), the function calculates the distribution factor ($AF$) to assess the spread of the obtained solution set. Meanwhile, Equation (35) defines the spatial performance indicator ($SP$), which evaluates the uniformity of the distribution of the solution set. A smaller $AF$ value suggests a broader spread of the solution set, whereas a smaller $SP$ value indicates a more uniform distribution. In this section, the proposed convergence metrics $CdI$, $EdI$, and $EdN$ are evaluated. The $AF$ and $SP$ values for the solution sets obtained by the algorithm under 2-, 3-, and 4-dimensional optimization objectives are illustrated in Figure 5, Figure 6 and Figure 7.

$$\begin{array}{cc}\hfill & AF=\sqrt{{\displaystyle \frac{1}{\left|A\right|-1}}\sum _{i=1}^{\left|A\right|}{\left(d_i^{\ast }-\overline{d}\right)}^2}\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & SP=\sqrt{{\displaystyle \frac{1}{\left|A\right|-1}}\sum _{i=1}^{\left|A\right|}{\left(\overline{d}-d_i\right)}^2}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & d_i^{\ast }=\underset{j\in A\setminus j\ne i}{min}\sum _{k=1}^m\left|f_k^i-f_k^j\right|,\overline{d}=\left(\sum _{i=1}^{\left|A\right|}{d}_i^{\ast }\right)/\left|A\right|,i,j\in A\hfill \end{array}$$

(36)

where m represents the number of optimization objective functions, $d_i$ denotes the minimum distance from the i-th solution to the other solutions in set A, and $\overline{d}$ is the average of all $d_i$ values.

Figure 5, Figure 6 and Figure 7 demonstrate that the solution sets generated by the MC-ObOEA algorithm exhibited favorable $AF$ and $SP$ values across all convergence indicators. This suggests that the solution sets are both extensive and uniformly distributed, making all three indicators effective for solving Pareto front optimization problems. Further analysis shows that for two optimization objectives, the $AF$ and $SP$ values were relatively larger across all network scales and convergence indicators, likely due to the smaller number of frontier solutions leading to a more sparse distribution among individuals. As the number of optimization objectives increased to three and four, the number of frontier solutions rose, resulting in decreased $AF$ and $SP$ values and, consequently, a more extensive and uniform distribution of the solution sets. For example, with four optimization objectives and a 15-node network, the $EdN$ indicator yielded the smallest $AF$ and $SP$ values, slightly outperforming the $CdI$ indicator. However, the difference between $EdN$ and $CdI$ was minimal—only 0.83—whereas the difference between the $EdI$ indicator and both $AF$ and $SP$ was significantly larger at 29. For the 30-node and 50-node networks, the $EdI$ indicator produced the smallest $AF$ and $SP$ values, indicating higher-quality solution sets, with only minor differences between the $CdI$ and $EdN$ indicators. In short, the MC-ObOEA algorithm performed effectively across all three convergence indicators, demonstrating robust results in optimizing solution set quality.

The solution set coverage indicator ($C\_metric$) assesses the dominance relationship between two solution sets, as mathematically defined in Equation (37).

Table 7. Comparison of solution sets $C\_metric$ under different convergence indicators.

Optimization Problem	$\mathit{CdI}$	$\mathit{EdI}$	$\mathit{EdI}$	$\mathit{EdN}$	$\mathit{CdI}$	$\mathit{EdN}$
15-2	0.0125	0.1250	0.1429	0.1250	0.1429	0.1250
15-3	0.3333	0.0345	0.0345	0.3333	0.0345	0.0345
15-4	0.0345	0.0345	0.0333	0.0345	0.0333	0.0345
30-2	0.0769	0.0769	0.0769	0.0769	0.0769	0.0769
30-3	0.0417	0.0345	0.0345	0.0417	0.0345	0.0345
30-4	0.0345	0.0385	0.0385	0.0345	0.0385	0.0385
50-2	0.2000	0.2000	0.1111	0.2000	0.1111	0.2000
50-3	0.3230	0.0400	0.0370	0.0323	0.0370	0.0400
50-4	0.0417	0.0278	0.0000	0.0417	0.0000	0.0278

Note: The “Optimization Problem” column denotes the format “number of nodes in the network—number of optimization objectives”. For instance, a problem labeled as 15-2 indicates an optimization problem involving a network with 15 nodes and 2 objectives.

presents the coverage rates for the solution sets obtained using the $CdI$, $EdI$, and $EdN$ indicators across different optimization problems.

$$C\_metric(A,B)={\displaystyle \frac{\left\{\left|\left\{u\in \left.B\right|\exists v\in A:v\mathrm{dominates}u\right\}\right|\right\}}{\left|B\right|}},$$

(37)

where the numerator represents the number of solutions in set B that are dominated by at least one solution from set A, while the denominator represents the total number of solutions in set B.

The data in Table 7 clearly shows that the solution sets derived from the $CdI$ and $EdI$ indicators exhibited minimal overlap across all optimization problems. A comparison between the $EdI$ and $EdN$ indicators indicates that $EdN$ only outperformed $EdI$ in the 50-node, 4-objective optimization problem. Similarly, when comparing $CdI$ and $EdN$, $EdN$ only surpassed $CdI$ in the same 50-node, 4-objective scenario. Overall, Table 7 highlights the effectiveness of all three convergence indicators in solving this test problem.

4.3. Impact of Distribution Indicators on the Performance of MC-ObOEA

We utilized distribution indicators based on cosine similarity ($Cs$) and Euclidean distance ($Ed$), running the algorithm 20 times to address optimization problems with two, three, and four objectives across 15-node, 30-node, and 50-node networks.

Table 8. Comparison of solution sets $AF$, $SP$, and $C\_metric$ under different methods.

Optimization Problem	$\mathit{AF}$		$\mathit{SP}$		$\mathit{C}\_\mathit{metric}$
	$\mathit{Cs}$	$\mathit{Ed}$	$\mathit{Cs}$	$\mathit{Ed}$	$\mathit{Cs}$	$\mathit{Ed}$
15-2	1272.79	1305.11	1266.09	1346.38	0.1670	0.1250
15-3	397.49	413.82	362.42	381.76	0.0340	0.0330
15-4	394.80	461.91	401.38	423.56	0.0303	0.0345
30-2	822.69	941.07	809.31	935.13	0.0769	0.0769
30-3	346.42	698.58	282.68	624.86	0.0476	0.0417
30-4	429.68	668.30	379.47	637.12	0.0333	0.0345
50-2	760.31	793.15	661.75	681.66	0.1250	0.2000
50-3	554.10	819.72	489.89	718.31	0.0476	0.0323
50-4	364.63	649.91	298.34	508.81	0.0556	0.0417

provides the mean values of the $AF$, $SP$, and $C\_metric$ indicators for the solution sets obtained through different methods. Figure 8 illustrates the variation in $AF$ and $SP$ values of the solution sets across evolutionary generations for the 15-node network.

From Table 8 and Figure 8, it is evident that the $Cs$-based method consistently produces smaller values $AF$ and $SP$ compared to the $Ed$-based method in the nine optimization problems. Although there was no significant difference in the values of $C\_metric$ between the two methods, the overall results suggest that the algorithm utilizing cosine similarity offers superior performance in terms of distribution metrics. Further analysis shows that as the number of objectives increased, the $Cs$ method outperformed the $Ed$ method in both the breadth and uniformity of the solution set distribution. For example, the $AF$ values for the $Cs$ method were lower than those for the $Ed$ method, with differences of 32.32, 16.33, 67.11; 118.38, 352.16, 238.62; and 32.84, 265.62, 285.28, respectively. The $SP$ values followed a similar trend. As illustrated in Figure 8, the $Cs$-based method reduced the $AF$ and $SP$ values more quickly in the early iterations and stabilized sooner. In conclusion, the $Cs$ distribution metrics provide more desirable performance for frontier solution sets than the Euclidean distance-based approach in this test problem.

To more intuitively demonstrate the effectiveness of using the $Cs$ method in the algorithm, we applied both methods to optimize the objectives for a 15-node network with four objectives, which are visualized using a parallel coordinate system. Figure 9 shows the solution sets obtained from a single run of the MC-ObOEA algorithm under each method, where each curve represents a different solution, and the points along the curves correspond to the objective values. To standardize the four objective values within the same coordinate system, data processing was performed on the objective values of each solution, with the original values of Objectives 1 to 4 scaled by factors of 6 × ${10}^4$, $1\times {10}^4$, $1\times {10}^3$, and 1, respectively. From Figure 9, it can be observed that the objective value distribution for the solution set obtained using the $Cs$ method was more uniform, whereas the $Ed$ method yielded a higher number of similar solutions, indicated by the overlapping curves. This shows that the $Cs$ method more effectively differentiates similar solutions, thus improving the overall quality of the solution set.

4.4. Algorithm Effectiveness Analysis

Based on the 15-node network established in this study and the specified algorithm parameters, Figure 10 illustrates the iterative evolution of the MC-ObOEA high-dimensional multiobjective algorithm. It depicts the average value of each objective for all individuals in the population across generations during the iteration process. As shown in Figure 10, the algorithm consistently converged toward minimizing individual objective values, confirming the correctness of the algorithm’s design.

Figure 11 shows the distribution of the objective values for a set of Pareto frontier solutions obtained by running the algorithm. The four box plots, from left to right, represent the total transportation time (in minutes), the total transportation cost (in CNY), the carbon emission cost (in CNY), and the rate of cargo loss. As depicted, the values of all four objective functions stabilized after a certain number of iterations, demonstrating the algorithm’s convergence and initially validating the effectiveness of its solution approach.

The algorithm was executed 20 times for each of the three different network scales, resulting in 20 sets of Pareto frontier solutions.

Table 9. Pareto frontier solutions for different network scales.

Network Scale	Stats	Obj.1	Obj.2	Obj.3	Obj.4
15-node	Min	6721.7	5184.4	301.5	0.21
	Mean	9637.52	5724.16	397.53	0.28
	Std	246.49	72.13	2.21	0.0072
30-node	Min	9633.35	7158.5	413.6	0.28
	Mean	14,483.72	7835.96	615.73	0.38
	Std	98.05	60.21	2.01	0.0024
50-node	Min	13,336	9694.5	609.5	0.29
	Mean	18,497.79	10,359.2	770.52	0.39
	Std	556.16	150.37	23.62	0.0068

presents the average minimum values and the overall mean values for each objective in these sets of solutions. For example, in the 15-node network, the average minimum values for total transportation time (Obj.1), total transportation cost (Obj.2), carbon emission cost (Obj.3), and cargo loss rate (Obj.4) are 6721.7, 5184.4, 301.5, and 0.21, respectively, with corresponding mean values of 9637.52, 5724.16, 397.53, and 0.28. The standard deviation, which reflects the variability in the minimum values of these objective functions across runs, is 246.49, 72.13, 2.21, and 0.0072. This indicates that the algorithm consistently achieves stable optimization performance across networks of varying scales.

To further validate the effectiveness of the proposed algorithm, the MC-ObOEA algorithm was compared with the RPEA [28] and NSGA-II [29] algorithms.

Table 10. Quality comparison of solution sets obtained by different algorithms with different evaluation indicators.

Evaluation Indicators	$\mathbf{AN}$			$\mathbf{AF}$			$\mathbf{SP}$
Optimization Problem	MC-ObOEA	RPEA	NSGA-II	MC-ObOEA	RPEA	NSGA-II	MC-ObOEA	RPEA	NSGA-II
15-2	6.3	7.6	6.4	1272.79	1072.72	1254.86	1266.09	1006.77	1250.99
15-3	7.6	11.3	14.25	397.49	418.41	473.14	362.42	398.51	438.57
15-4	9.3	14.4	17.8	393.8	428.69	455.13	405.38	393.26	419.62
30-2	8.1	8.5	9.15	826.69	638.64	815.04	809.31	565.41	788.012
30-3	11.6	14.3	16.55	346.42	476.95	603.94	282.68	419.25	567.8
30-4	12.9	15.5	19	427.68	540.92	566.52	379.47	382.17	526.98
50-2	8.7	9	10.5	754.3	722.86	742.28	661.74	693.01	699.97
50-3	11.6	14.5	19.7	554.1	534.46	524.25	489.89	496.1	490.45
50-4	13.4	16.3	23.4	364.63	410.13	493.6	298.34	433.19	415.18

presents the average number of frontier solutions ($AN$), the average distribution ($AF$), and the values of the uniformity indicator ($SP$) obtained by the different algorithms in nine optimization problems after 20 runs. The paired samples t-test analysis allows us to assess the degree of difference between the data by examining the size of the p-value within the confidence interval. If the p-value is less than 0.05, it indicates a significant difference; otherwise, no significant difference is present.

Table 11. Result of paired sample t-test for solution sets obtained by different algorithms.

Evaluation Indicator	Algorithm	t-Value	p-Value	$\mathit{\alpha }$ = 0.05
$AN$	MC-ObOEA-PREA	4.710	0.002	yes
	MC-ObOEA-NSGA-II	4.453	0.002	yes
	RPEA-NSGA-II	3.459	0.009	yes
$AF$	MC-ObOEA-PREA	0.268	0.795	no
	MC-ObOEA-NSGA-II	2.042	0.075	no
	RPEA-NSGA-II	3.226	0.012	yes
$SP$	MC-ObOEA-PREA	0.392	0.706	no
	MC-ObOEA-NSGA-II	2.154	0.063	no
	RPEA-NSGA-II	2.666	0.029	yes

presents the results of the paired samples t-test for different algorithms, which were calculated using SPSS based on the data from Table 10. Additionally,

Table 12. Quality comparison of solution sets obtained by different algorithms.

Optimization Problem	MC-ObOEA	RPEA	MC-ObOEA	NSGA-II	RPEA	NSGA-II
15-2	0.1429	0.1111	0.1250	0.1111	0.1670	0.1000
15-3	0.0357	0.0357	0.0313	0.0357	0.0370	0.0300
15-4	0.0357	0.0294	0.0323	0.0000	0.0370	0.0290
30-2	0.0769	0.0625	0.0769	0.0667	0.1000	0.0380
30-3	0.0400	0.0385	0.0500	0.0385	0.0500	0.0300
30-4	0.0345	0.0370	0.0345	0.0333	0.0430	0.0320
50-2	0.0270	0.0286	0.1667	0.1250	0.1670	0.1100
50-3	0.0345	0.0370	0.0345	0.0333	0.0450	0.0260
50-4	0.0323	0.0000	0.0417	0.0333	0.0360	0.0230

lists the mutual coverage rates of the solution sets produced by each algorithm. From Table 12, it is clear that the mutual dominance between the solution sets of the three algorithms was rare, indicating no significant differences in this regard. Although Table 11 indicates that the $AF$ and $SP$ values of MC-ObOEA for the five optimization problems are not significantly different from those obtained with PREA and NSGA-II, the results in Table 10 show that MC-ObOEA achieved the lowest $AF$ and $SP$ values in five out of the nine optimization problems, demonstrating its superior overall performance. Additionally, the AN value of the MC-ObOEA was tested against the other two algorithms, yielding a p-value of 0.002, which is less than 0.05. This indicates a statistically significant difference, highlighting MC-ObOEA’s effectiveness in enhancing the diversity of the solution set. The $AN$ values for MC-ObOEA ranged from 6.3 to 13.4, with a more gradual increase as the number of optimization objectives increased. In contrast, the NSGA-II algorithm exhibited a sharp rise in the number of frontier solutions as both network size and the number of optimization objectives increase, indicating a degree of diversity in its solution sets. However, as shown in Table 12, the solution sets obtained by the NSGA-II algorithm did not exhibit a significant dominance over those produced by the other two algorithms, suggesting that the NSGA-II algorithm generates a substantial number of similar solutions.

Using the example of the 30-node network with three optimization problems, when the number of objectives is two, the differences in $AN$ values between the MC-ObOEA algorithm and the RPEA and NSGA-II algorithms are −0.4 and −1.05, respectively, indicating no significant variation. However, the differences in $AF$ values are 188.05 and 11.65, while the differences in $SP$ values are 243.9 and 21.298, suggesting that the RPEA and NSGA-II algorithms achieved higher solution set quality and better performance in this scenario. When the number of objectives increased to three and four, significant differences emerged. The $AN$ value differences between MC-ObOEA and the other two algorithms became −2.7, −4.95 (for three objectives) and −2.6, −6.1 (for four objectives). The $AF$ value differences are −130.53, −257.52 (for three objectives) and −113.24, −138.84 (for four objectives). Similarly, the $SP$ value differences were −136.57, −285.12 (for three objectives) and −2.7, −147.51 (for four objectives), all of which highlight substantial discrepancies. These results indicate that the MC-ObOEA algorithm outperformed both the RPEA and NSGA-II algorithms in handling high-dimensional multiobjective optimization problems.

4.5. Effects of Parameter Variations on Carbon Emissions and Cargo Loss

The carbon tax rate is a key factor in determining carbon emission costs, as fluctuations in the tax rate can significantly affect the algorithm’s optimization process. To assess this impact, a 15-node network was used as an example, with four carbon tax rates tested: 10 CNY$/t$, 25 CNY$/t$, 40 CNY$/t$, and 60 CNY$/t$. The algorithm was run for each tax rate, and transportation plans with identical routes and modes were selected from the Pareto solution set. These were then compared to analyze the changes in carbon emission costs under different tax rates. The results are presented in

Table 13. Transport options under different carbon tax rates.

Routes and Modes	Carbon Tax Rate (CNY/t)	Carbon Emissions Cost (CNY)	Rate of Change	Routes and Modes	Carbon Tax Rate (CNY/t)	Carbon Emissions Cost (CNY)	Rate of Change
Scheme 1: o(H)C(H) D(H)K(H)d	10	375	/	Scheme 2: o(W)C(R) D(R)K(R)d	10	274	/
	25	520	0.37		25	300	0.095
	40	660	0.27		40	326	0.087
	60	850	0.29		60	360	0.104
Scheme 3: o(H)A(H) E(R)M(R)d	10	326	/	Scheme 4: o(H)A(R) E(R)M(R)d	10	300	/
	25	401	0.23		25	357	0.19
	40	476	0.19		40	405	0.13
	60	578	0.21		60	469	0.16
Scheme 5: o(R)B(H) D(H)K(R)d	10	340	/	Scheme 6: o(R)B(W) D(R)K(R)d	10	272	/
	25	422	0.24		25	298	0.095
	40	504	0.19		40	323	0.084
	60	611	0.21		60	357	0.11
Scheme 7: o(H)B(H) G(H)L(R)d	10	352	/	Scheme 8: o(R)B(R) G(H)L(R)d	10	310	/
	25	455	0.29		25	358	0.15
	40	558	0.23		40	407	0.14
	60	695	0.25		60	472	0.16

.

Table 13 presents eight transportation schemes, including all-road transportation, rail–water intermodal transportation, and various mixed modes. Comparing schemes 1 and 2, it is evident that the rate of increase in carbon emission costs for scheme 1 was consistently higher than for scheme 2 as the carbon tax rate rose. For instance, when the carbon tax rate increased from 10 CNY$/t$ to 60 CNY$/t$, the carbon emission cost for scheme 1 increased by 475 CNY, while for scheme 2, the increase was only 86 CNY. This highlights the importance of reducing long-distance road transportation to reduce carbon emissions and move towards greener modes such as rail and waterway transport. A similar trend was seen when comparing schemes 7 and 8. Scheme 7, which relies more heavily on road transportation, showed greater fluctuations in carbon emission costs compared to scheme 8. However, comparing schemes 2 and 6, both of which utilize rail–water intermodal transportation, reveals that carbon emission costs were less sensitive to changes in the carbon tax rate, suggesting that such schemes are more resilient to tax fluctuations. These comparisons demonstrate that increasing the carbon tax rate leads to higher carbon emission costs across transportation schemes. Consequently, freight companies are incentivized to favor cleaner, greener modes of transportation, such as railways and waterways.

The cold chain and fresh goods are highly perishable and sensitive to spoilage, requiring strict control over the transport environment, particularly with respect to temperature fluctuations. The maintenance of stable temperatures during transport is vital for minimizing product loss. To better understand the risks of spoilage, it is essential to analyze the loss rates of various goods under different temperature conditions. Different types of products have varying activation energies, which represent the minimum energy required to initiate a chemical reaction, typically leading to spoilage. This activation energy, denoted $E\alpha$ and measured in kJ/mol, varies between products. Goods with lower activation energies are more prone to spoilage and degradation.

Table 14. Product activation energy table.

Product	Activation Energy (kJ/mol)	Product	Activation Energy (kJ/mol)	Product	Activation Energy (kJ/mol)
Apricot	30.62	Date	54.51	Tomato	32.94
Sliced potato	39.49	Bean	27.71	Purple grape	67.29
Sliced pumpkin	78.93	Garlic	23.48	Sliced carrot	25.93

provides a list of common items along with their respective activation energies.

Using the 15-node network as an example, the activation energy ($E\alpha$) was set at 34 kJ/mol, the temperature fluctuation ($\Delta F$) during transitions was set at 1 °C, the gas constant (R) was 8.314, and the maximum reaction rate ($K_{max}$) was $5\times {10}^4$. All other algorithmic parameters remained unchanged. Four transportation schemes were selected from the Pareto frontier solutions.

• Scheme 1: $o\left(H\right)B\left(H\right)G\left(H\right)L\left(R\right)d$, with a total transportation time of 7126 min, a transportation cost of 5857 CNY, a carbon cost of 490 CNY, and a cargo depletion rate of 0.220.
• Scheme 2: $o\left(R\right)A\left(H\right)E\left(R\right)M\left(R\right)d$, with a total transportation time of 8720 min, a transportation cost of 6033 CNY, a carbon cost of 375 CNY, and a cargo depletion rate of 0.263.
• Scheme 3: $o\left(W\right)C\left(R\right)D\left(H\right)K\left(W\right)d$, with a total transportation time of 15,565 min, a transportation cost of 5046 CNY, a carbon cost of 388 CNY, and a cargo depletion rate of 0.426.
• Scheme 4: $o\left(R\right)B\left(W\right)D\left(R\right)K\left(R\right)d$, with a total transportation time of 13,002 min, a transportation cost of 5484 CNY, a carbon cost of 306 CNY, and a cargo loss rate of 0.368.

The temperature variation ($\Delta F$) ranges from 0 to 2 °C, while the activation energy ($E\alpha$) varies between 30 and 48 kJ/mol.

Table 15. Cargo loss changes with activation energy and temperature difference for Scheme 1.

	$\mathit{E}\mathit{\alpha }$	30 kJ/mol	33 kJ/mol	36 kJ/mol	39 kJ/mol	42 kJ/mol	45 kJ/mol	48 kJ/mol
ΔF (°C)
0.2		0.72071	0.31413	0.08855	0.02513	0.00779	0.00212	0.00058
0.4		0.72522	0.31459	0.08917	0.02602	0.00775	0.00777	0.00058
0.6		0.73096	0.31562	0.09312	0.02624	0.00779	0.00776	0.00058
0.8		0.7355	0.31685	0.09641	0.02709	0.00779	0.00781	0.00058
1		0.7409	0.31727	0.09716	0.02816	0.00782	0.00777	0.00059
1.2		0.74419	0.31738	0.09942	0.02874	0.00781	0.00782	0.00059
1.4		0.75269	0.31828	0.09978	0.02884	0.00782	0.0078	0.00058
1.6		0.75361	0.31682	0.09925	0.02895	0.0078	0.0078	0.00059
1.8		0.75399	0.31752	0.09941	0.02899	0.00782	0.00787	0.00058
2		0.75436	0.31855	0.10007	0.029	0.00785	0.00777	0.00058

illustrates the cargo loss rate for transport Scheme 1 across different activation energy and temperature fluctuation levels. Figure 12a–d show the trends of variation in the cargo loss rates for the first to fourth transport schemes, respectively. In Figure 12, the activation energy ($E\alpha$) axis is incremented by 3 kJ/mol, while the temperature fluctuation axis ($\Delta F$) is increased by 0.2 °C.

Figure 12 reveals a general trend in which the cargo loss rate decreases as the activation energy increases, but increases with increased temperature fluctuations. As shown in Figure 12a and Table 15, under the same temperature conditions, the cargo loss rate was particularly sensitive to changes in activation energy. For example, when the activation energy reached 48 kJ/mol, the cargo loss rate dropped to 0.00058. Conversely, for the same cargo, an increase in temperature difference resulted in a higher loss rate, especially for perishable goods with lower activation energy. The impact of temperature fluctuations was more pronounced on goods with lower activation energy, underscoring the critical importance of controlling temperature during transport. For goods with low activation energy, such as fresh or perishable items with $E\alpha$ equal to 30 kJ/mol, the cargo loss rate remained elevated across the four transport schemes, particularly in Schemes 3 and 4. This highlights the need for the development of faster transport modes, such as air transport, to reduce the transportation time for sensitive goods in the future.

5. A Case Study

Assume a standard container shipment is to be transported from Chongqing, China, to Duisburg, Germany. The carrier must develop a feasible transportation plan that meets the requirements of the carrier while ensuring timely and high-quality delivery to the consignee. To achieve this, a transportation network is constructed with Chongqing as the starting point and Duisburg as the destination. The network incorporates various modes of transport, including road, rail, and water, and spans 24 city nodes, such as Chongqing, Xinjiang Alashankou, Shanghai, and Duisburg, as illustrated in Figure 13.

In Figure 13, nodes 1 to 24 represent cities or countries along the route: Chongqing, Xinjiang Alashankou, Shanghai, Ningbo-Zhoushan, Guangzhou, Beibu Gulf (Guangxi), Hong Kong, the Philippines, Singapore, Malaysia, Vietnam, Kazakhstan, Turkmenistan, Russia, Iran, India, Sri Lanka, Djibouti, Egypt, Belarus, Poland, France, Italy, and Duisburg, Germany. The specific port cities are determined by the actual logistics conditions. The transport modes are indicated as follows: H for the road, R for rail, and W for water. Lines connecting the nodes indicate their accessibility. The transportation times between nodes follow a normal distribution, with mean and standard deviation values presented in

Table 16. Multimodal transportation network data.

1–2;∼;[67,8]11,055;∼;	9–18;∼;∼;[342,32]14,179;
1–3;[96,5]8678;∼;[240,24]2972;	10–9;∼;∼;[28,8]628;
1–4;[80,5]6720;[62,6]4520;∼	10–16;∼;∼;[97,8]4097;
1–5;[89,7]8640;[60,7]5940;∼;	10–17;∼;∼;[76,8]3013;
1–6;∼;[14,3]4150;∼;	10–18;∼;∼;[196,18]9606;
2–12;∼;[82,5]10,276;∼;	11–9;∼;∼;[146,12]2246;
3–7;∼;∼;[84,24]3894;	11–10;∼;∼;[185,17]2585;
3–8;∼;∼;[106,8]2880;	12–13;∼;[68,4]2680;∼;
3–9;∼;∼;[192,24]5723;	12–14;∼;[32,2]2981;∼;
4–7;∼;∼;[37,8]1550;	13–15;∼;[32,2]2369;∼;
4–9;∼;∼;[127,21]4050;	14–20;∼;[12,4]3390;∼;
5–7;∼;∼;[96,8]2149;	15–18;∼;∼;[12,4]3689;
5–9;∼;∼;[125,14]5800;	16–18;∼;∼;[156,21]5568;
6–7;∼;∼;[26,8]3030;	17–18;∼;∼;[196,18]6196;
6–9;∼;∼;[126,18]6260;	18–19;∼;∼;[96,8]3109;
6–11;∼;∼;[196,22]6800;	19–22;∼;∼;[146,19]4146;
7–9;∼;∼;[129,8]3889;	19–23;∼;∼;[96,8]3533;
8–9;∼;∼;[97,18]3534;	20–21;∼;[12,2]1600;∼;
8–10;∼;∼;[123,26]3567;	21–24;∼;[42,6]6342;∼;
9–16;∼;∼;[172,12]4572;	22–24;[26,6]3916;[32,3]2890;[52,7]1820;
9–17;∼;∼;[136,16]3581;	23–24;[21,5]4211;[35,3]3501;[64,17]2671;

. The schedule of transport modes between nodes is homogeneous, and transshipment and reloading data remain unchanged at each node. The algorithm and model parameters are applied as previously described, with the solutions of the MC-ObOEA algorithm for this multimodal transportation network provided in

Table 17. Pareto frontier solution set obtained by the MC-ObOEA algorithm.

Scheme No.	Route and Mode	Obj.1	Obj.2	Obj.3	Obj.4
1	1(R)2(R)12(R)14(R)20(R)21(R)24	16,345	39,560	498	0.215
2	1(R)6(W)9(W)16(W)18(W)19(W)23(R)24	45,277	37,344	697	0.496
3	1(H)4(W)9(W)16(W)18(W)19(W)23(H)24	48,270	38,594	957	0.519
4	1(R)6(W)9(W)16(W)18(W)19(W)23(W)24	49,315	36,525	716	0.526
5	1(R)4(W)9(W)17(W)18(W)19(W)23(W)24	52,634	34,504	730	0.549
6	1(W)3(W)9(W)18(W)19(W)23(R)24	63,297	38,813	585	0.612

.

Table 17 presents six transport schemes included in the set of Pareto frontier solutions. Scheme 1 follows the Chongqing–Xinjiang–Europe route, which utilizes rail transport exclusively. It starts at Chongqing Tuanjiecun Rail Container Center Station and passes through Xinjiang Alashankou, Kazakhstan, Russia, Belarus, and Poland before reaching Duisburg, Germany. Schemes 2 and 4 share similar main nodes. Both start at the Chongqing Tuanjiecun Rail Container Center Station and use rail transport through the New Western Land-Sea Corridor to the Beibu Gulf Port in Nanning, Guangxi, from where they are shipped by water to Singapore, India, Djibouti, Egypt, and Italy. In Scheme 2, the journey continues by rail from Italy to Duisburg, while in Scheme 4, the remaining portion is completed by water. In Scheme 3, the domestic segment is transported by road from Chongqing to Ningbo-Zhoushan Port, followed by water for the international segment, with the final leg from Italy to Duisburg completed by rail. Scheme 5 employs the newly opened Yuyong channel of the river railway (since 2019) to transport cargo domestically to Ningbo Zhoushan Port, and the international section is completed primarily by waterways. Finally, Scheme 6 uses Chongqing’s traditional route to the sea, starting from Chongqing Guoyuan Port along the Yangtze River Golden Waterway to Shanghai Yangshan International Container Port. The international leg of this route then passes through Singapore, Djibouti, Egypt, and Italy before reaching Duisburg. Figure 14 illustrates the routes for each transport scheme.

Further analysis reveals that the primary advantage of Transport Scheme 1 is its relatively short transportation time of approximately 11 days, which minimizes cargo spoilage. Additionally, rail transport is an environmentally friendly transport mode, resulting in lower carbon emissions and reduced carbon costs. However, Scheme 1 has the highest total transportation cost compared to the other options. The cost of the Yuxin-Europe route could be optimized by streamlining customs clearance and minimizing the frequent transitions and tracking changes that occur during transport. Enhancing the efficiency of customs procedures and reducing intermediate transition costs offer potential areas for improvement.

Transport Scheme 5 has the lowest overall transportation cost. This is largely due to the use of the Chongqing-Ningbo high-speed rail corridor, which reduces domestic transport distance. In addition, preferential policies promoting rail-water intermodal transport help to lower costs. The international leg, which relies entirely on water transport, contributes further to this cost efficiency. In contrast, Transport Scheme 3 incorporates road transport, leading to the highest carbon emissions cost of 957 CNY, highlighting its significant environmental impact.

Transport Schemes 2 and 4 show intermediate performance in all four objective functions. Scheme 2 uses rail transport from Italy to Duisburg, while Scheme 4 opts for water transport for the same segment. Although Scheme 2 offers shorter transit times, it incurs higher transportation costs, suggesting that water transport is more cost-effective, whereas rail transport provides greater time efficiency.

Transport Scheme 6, though characterized by relatively low carbon emissions, does not perform notably well in other areas. Its domestic segment involves water transport, resulting in longer transit times and higher cargo loss rates. Additionally, its total transportation cost remains relatively high, making it less favorable for perishable goods. However, for bulk goods that are not sensitive to time, Scheme 6 remains a viable option due to its established and reliable logistics operations.

6. Conclusions

This study developed a high-dimensional, multiobjective route optimization model for multimodal cold chain transportation, addressing key objectives such as minimizing total transportation costs, transit time, carbon emissions, and food loss rates. The model integrates a comprehensive range of factors, including en-route costs, trans-shipment and storage costs, and uncertainties in transportation time, using Monte Carlo simulation to enhance robustness. A carbon pricing function quantifies carbon emissions, while a Weibull function models the loss rate of perishable goods as a function of temperature and time, reflecting the real-world complexities of cold chain logistics.

To solve this intricate optimization problem, we proposed the MC-ObOEA, a multi-objective evolutionary algorithm that combines Monte Carlo simulation with a one-by-one selection strategy. The MC-ObOEA effectively balances convergence and diversity indicators, producing high-quality Pareto frontier solutions. Evaluation using three convergence metrics ($CdI$, $EdI$, $EdN$) demonstrates the stability and efficiency of the algorithm to solve multidimensional high-dimensional optimization problems. In particular, the $Cs$ method outperformed the $Ed$ method in terms of solution distribution and uniformity, underscoring the algorithm’s capability to generate diverse and well-distributed solutions across the Pareto frontier. The MC-ObOEA consistently exhibited strong optimization performance, particularly in managing the complexities of high-dimensional, multi-objective problems. The parameter impact analysis further highlights the sensitivity of carbon emission costs to carbon tax rates and the critical role temperature variation plays in the loss rates of perishable goods, especially those with low activation energy. These insights underscore the importance of precise temperature control in cold chain logistics to minimize food spoilage and improve transportation efficiency.

This study has several limitations. First, the parameter values used were adopted directly from previous literature without a detailed analysis of their selection process or an examination of how different parameter configurations might impact algorithm performance and the solution set. Future research could include sensitivity analysis, stress testing, and robustness analysis to strengthen the justification for these parameter choices. Additionally, comparative analysis with more advanced algorithms is limited, suggesting that future work could explore integrating and comparing our approach with other cutting-edge algorithms to further enhance performance. Although this study focuses on optimizing single-source to single-destination multimodal transportation routes, future research should expand the model to accommodate multisource to multidestination scenarios. This extension offers promising avenues for further exploration in the optimization of cold chain logistics, with the potential to improve efficiency and sustainability in global food distribution networks.