A Multi-Objective Model for Sustainable Perishable Food Distribution Considering the Impact of Temperature on Vehicle Emissions and Product Shelf Life

: The food distribution process is responsible for signiﬁcant quality loss in perishable products. However, preserving quality is costly and consumes a tremendous amount of energy. To tackle the challenge of minimizing transportation costs and CO 2 emissions while also maximizing product freshness, a novel multi-objective model is proposed. The model integrates a vehicle routing problem with temperature, shelf life, and energy consumption prediction models, thereby enhancing its accuracy. Non-dominated sorting genetic algorithm II is adapted to solve the proposed model for the set of Solomon test data. The conﬂicting nature of these objectives and the sensitivity of the model to shelf life and shipping container temperature settings are analyzed. The results show that optimizing freshness objective degrade the cost and the emission objectives, and the distribution of perishable foods are sensible to the shelf life of the perishable foods and temperature settings inside the container.


Introduction
Improving the sustainability of a distribution network involves tradeoffs between multiple conflicting objectives, including minimizing transportation costs (e.g., fuel and vehicle maintenance costs, driver salaries), fulfilling customer requirements (e.g., on-time deliveries, short lead times), and limiting environmental impact (e.g., vehicle emissions). However, optimizing sustainability in perishable food distribution is particularly challenging, primarily because of temperature control requirements [1]. Temperature is a major determinant of the shelf life of a perishable product [2][3][4]. Even small or infrequent deviations from recommended temperature settings can significantly reduce product shelf life [5][6][7] because increased temperature accelerates the growth rate of the microorganisms that are responsible for quality degradation in perishable foods [5,8]. Although refrigerated vehicles' cargo is well-isolated, it can experience frequent exposure to increased temperature when the vehicle stops to make deliveries to other customers [2,9]. As a result, an estimated 8-23% loss in perishable food quality occurs during the distribution process [10].
This loss in quality increases the likelihood that the food is wasted [11]. According to the United States Department of Agriculture, 30-40% of food in the U.S. is wasted [12], with 40% of these losses occurring post-harvest [13]. As a result, much of the resources consumed by the production and distribution of perishable food, as well as their associated environmental impacts, are in vain [13]. It is estimated that food waste costs the U.S. economy $218 billion each year [14]. Moreover, as described

Literature Review
There is a rich literature related to the distribution of perishable food. The review presented in this paper focuses on literature that uses mathematical modeling to optimize food distribution systems.
Transit time and temperature are the two most influential factors on the quality of delivered perishable food products. Integrating product transit time is straightforward for mathematical models that are already tracking delivery times for other purposes, such as time window constraints. Therefore, many studies that take perishability into account consider the reduction of transit time directly or indirectly in their models [10,[26][27][28][29][30][31]. These models assume that the rate of product deterioration is constant over time, such that minimizing product time in transit is linearly equivalent to an increase in the quality of the delivered products [32]. Hence, transit time is used as a proxy for product quality loss, either as part of cost minimization or revenue maximization objective.
However, the inherent tradeoffs between delivering high-quality products and minimizing distribution costs have inspired researchers to consider them as separate objectives. The multi-objective model developed by Bortolini et al. [32] minimizes delivery time as a separate objective to represent the freshness of delivered foods. Similarly, Amorim and Almada-Lobo [33] Musavi and Bozorgi-Amiri [25], and Rahbari et al. [23] used a multi-objective approach to study product freshness maximization, where freshness was linearly estimated. Amorim, Günther, & Almada-Lobo [34] and Hsu, Chen, and Wu [35] integrated the effect of temperature into the rate of quality degradation, such that the rate of quality degradation increases with increasing storage temperature. Since the growth rate of microorganisms that spoil food is exponential, the impact of temperature variations on product shelf life can be estimated using exponential functions, which can then be maximized [22,[36][37][38].
Temperature-controlled transportation is energy-intensive and consequently releases a large volume of CO 2 emissions into the environment. The energy required to carry a load between two locations is evaluated in Bektaş & Laporte's [39] study with respect to traveled distance, vehicle acceleration speed, the slope of the road, the weight of the load, the density of air, and the frontal surface of the vehicle. Stellingwerf et al. [17] adapted this methodology to evaluate CO 2 emissions for both transportation and refrigeration, based on the assumption of constant energy losses during unloading and through the walls of the vehicle. Wang et al. [38] transformed CO 2 emissions and the energy required for transportation and refrigeration into costs, which were minimized. Accorsi, Gallo, & Manzini [40] calculated the energy consumption of distribution activities, which were then minimized in the objective function. Hsu et al. [31] integrated the effects of energy required for refrigeration of perishable products as a part of a cost minimization objective, where the cost of energy is a function of the constant and predetermined temperature inside the container, ambient temperature, volume of the container, duration of each stop, and frequency of opening the container. Hsu et al. [35] optimized the delivery of perishable products having different temperature requirements, accounting for refrigeration costs as a component of distribution costs. Some studies include CO 2 emissions in a separate objective of a multi-objective model (see e.g., Bortolini et al. [32], Govindan [25], and Wang et al. [43]), allowing decision-makers to assess the impact of reducing CO 2 emissions on other food distribution system objectives, as well as the marginal cost of reducing environmental effects.
Most models in the literature seek to maximize product quality throughout the distribution process. By contrast, Devapriya et al. [44], Khalili-Damghani et al. [45] and Nakandala, Lau, & Zhang [46] added a set of constraints to their model to ensure that the quality of the delivered products meets customers' expectations. Hsu et al. [31] calculated the volume of spoiled products based on the duration and ambient temperature of each delivery stop. Novaes et al. [9] used commercial software to predict the temperature inside the shipping container, using time-temperature data to evaluate the quality of products at each delivery location via a statistical indicator in a traveling salesman problem.
The impact of sustainability goals was reflected in multi-objective models. Several studies addressed the environmental and economic effects of the food supply chain in bi-objective models (see e.g., Govindan et al. [41], Soysal et al. [47], and Validi et al. [48]). The quality and safety of Sustainability 2020, 12, 6668 4 of 21 perishable products are affected by storage and distribution process time and temperature. Hence, Accorsi, Baruffaldi, and Manzini [49] evaluated the impact of temperature on operation efficiency and the safety of perishable products in a bi-objective study. Bortolini et al. [32] added the delivery time minimization as the third objective to the environmental and economic objectives to reduce the effect of traveling time on losing the quality of perishable products.
A diverse set of solution strategies have been used to solve the perishable food distribution problem. Ghezavati et al. [27] adapted a Benders decomposition model to solve a mixed-integer linear program. Chen et al. [26] and Farahani et al. [30] adapted a heuristic algorithm to solve the distribution planning problem as part of their production and distribution model. Metaheuristic approaches have also been widely used to find good-quality solutions in a reasonable amount of time. These approaches are also useful for finding a Pareto optimal frontier in multi-objective problems. Musavi and Bozorgi-Amiri [25] applied NSGA-II on a multi-objective hub location scheduling problem, in which total transportation costs and carbon emissions were minimized and food freshness was maximized. Amorim and Almada-Lobo [33] applied ε-constraint method to a small-scale problem and NSGA-II to a large-scale multi-objective problem that aimed to minimize total routing costs and maximize average freshness in a food distribution problem. Khalili-Damghani et al. [45] solved a bi-objective location-routing problem for the distribution of perishable products using the ε-constraint and NSGA-II algorithm. Their result showed that the NSGA-II solutions were as good as the ε-constraint solutions, but the metaheuristic algorithm outperformed the exact method in terms of solving time, especially for larger-scale problems. Govindan et al. [41] used a hybrid approach that integrated an adapted multi-objective particle swarm optimization (MOPSO) and an adapted multi-objective variable neighborhood search (AMOVNS) to solve a bi-objective location routing problem for perishable products with economic and environmental minimization objectives.
A summary of the important features of the related literature is given in Table 1. The related literature is ordered from oldest to the newest, and "x" represents that the publication included the specific feature in their study. Most of these studies either focus on food product perishability or the environmental impacts of temperature-controlled distribution. Only two studies cover cost, freshness, and emissions, and only Hsu et al. [31] and Novaes et al. [9] have considered the effect of temperature on product quality and carbon emissions.
To the best of our knowledge, the multi-objective vehicle routing problem (VRP) model presented in this paper is the first to use an integrated temperature prediction method to estimate product quality and refrigeration energy consumption while accounting for cost, freshness, and emissions in a perishable food distribution system. Specifically, this paper extends the MO-SVRP by adding a heat exchange model to accurately estimate the temperature inside the refrigerated container. This allows for a more accurate prediction of product freshness (i.e., shelf life) upon delivery, as well as improving the estimation of total emissions generated by refrigerated trucks. This paper also utilizes a novel adaptation of the NSGA-II metaheuristics algorithm to solve the MO-SVRP model.

Problem Statement and Model Formulation
Distribution problems typically seek to consolidate the flow of goods from a depot to their demand destinations into fewer routes [50]. The vehicle routing problem (VRP) is often utilized to formulate this problem, such that an optimal (i.e., minimum distance) route is determined, subject to constraints such as route connectivity, vehicle capacity limits, and the number of available vehicles. The model described in this paper integrates a VRP with methods that accurately predict temperature, estimate CO 2 emissions, and predict product freshness.
The proposed MO-SVRP model consists of the interconnected perishable food distribution modules including distribution model, temperature prediction, CO 2 emission, and shelf life prediction. The distribution model provides vehicles' delivery sequence and traveled distance which then are used to predict temperature inside the vehicle. The predicted temperatures are the key element of estimating the shelf life of the foods at their destination. The predicted temperature along with routing information enables the CO 2 estimation module to predict the environmental effect for the proposed routing plan. The predicted freshness of perishable products, CO 2 emission estimates, and distribution costs are the sustainability objectives of the MO-SVRP model. Figure 1 presents the integrated structure of the MO-SVRP model. In the remainder of this section, each of these methods and their mathematical relations are presented, and then the MO-SVRP assumptions, notations, and model are illustrated.

Problem Statement and Model Formulation
Distribution problems typically seek to consolidate the flow of goods from a depot to their demand destinations into fewer routes [50]. The vehicle routing problem (VRP) is often utilized to formulate this problem, such that an optimal (i.e., minimum distance) route is determined, subject to constraints such as route connectivity, vehicle capacity limits, and the number of available vehicles. The model described in this paper integrates a VRP with methods that accurately predict temperature, estimate CO2 emissions, and predict product freshness.
The proposed MO-SVRP model consists of the interconnected perishable food distribution modules including distribution model, temperature prediction, CO2 emission, and shelf life prediction. The distribution model provides vehicles' delivery sequence and traveled distance which then are used to predict temperature inside the vehicle. The predicted temperatures are the key element of estimating the shelf life of the foods at their destination. The predicted temperature along with routing information enables the CO2 estimation module to predict the environmental effect for the proposed routing plan. The predicted freshness of perishable products, CO2 emission estimates, and distribution costs are the sustainability objectives of the MO-SVRP model. Figure 1 presents the integrated structure of the MO-SVRP model. In the remainder of this section, each of these methods and their mathematical relations are presented, and then the MO-SVRP assumptions, notations, and model are illustrated.

Shelf life Prediction
Delivery Sequence

Temperature Prediction Based on Heat Transfer
The cooling unit in a refrigerated vehicle is constantly trying to preserve the temperature inside the container by blowing cold air. However, each time the vehicle makes a delivery stop and opens the container door for unloading, the heat exchange between the hot ambient air and cold container air raises the temperature inside the container. The capacity of the cooling equipment and the amount of heat exchange are the main factors that determine the temperature inside the container. In the MO-SVRP model, energy balance equations are applied to predict the temperature inside the container.
Since heat exchange when the vehicle is in transit differs from unloading, the container temperature in each of these stages is predicted using different methods. When the vehicle is moving, the cooling system runs until the temperature inside the container reaches the desired level (Td), and a thermostat turns the cooling system off. In most cases, the vehicle's cooling system and engine are turned off during unloading to protect the engine and to avoid polluting the air around the delivery dock. Therefore, it is assumed that the cooling system is not running during unloading, and therefore the temperature inside the container can only increase.

Temperature Prediction during Unloading
According to the energy balance equation, the overall heat that enters the container is equal to the accumulated heat inside the container (Equations (1)-(2); Holdsworth et al. [51]). The heat entering and accumulating inside the container at customer location j are denoted by AEj and AHj, respectively.

Temperature Prediction Based on Heat Transfer
The cooling unit in a refrigerated vehicle is constantly trying to preserve the temperature inside the container by blowing cold air. However, each time the vehicle makes a delivery stop and opens the container door for unloading, the heat exchange between the hot ambient air and cold container air raises the temperature inside the container. The capacity of the cooling equipment and the amount of heat exchange are the main factors that determine the temperature inside the container. In the MO-SVRP model, energy balance equations are applied to predict the temperature inside the container.
Since heat exchange when the vehicle is in transit differs from unloading, the container temperature in each of these stages is predicted using different methods. When the vehicle is moving, the cooling system runs until the temperature inside the container reaches the desired level (Td), and a thermostat turns the cooling system off. In most cases, the vehicle's cooling system and engine are turned off during unloading to protect the engine and to avoid polluting the air around the delivery dock. Therefore, it is assumed that the cooling system is not running during unloading, and therefore the temperature inside the container can only increase.

Temperature Prediction during Unloading
According to the energy balance equation, the overall heat that enters the container is equal to the accumulated heat inside the container (Equations (1)-(2); Holdsworth et al. [51]). The heat entering and accumulating inside the container at customer location j are denoted by AE j and AH j , respectively.
m a : Air mass (kg) r a : Air transfer ratio s a : Specific heat of the air (J kg −1 K −1 ) T j : The ambient temperature at location j (K) T 0 : Current temperature inside the container (K) t: Portion of unloading time (t ≤ u j ) m j : Cargo mass at location j (kg) m a : Air mass (kg) s c : Specific heat of cargo (J kg −1 K −1 ) s a : Specific heat of the air (J kg −1 K −1 ) dT dt : Rate of change in temperature (K) Following the energy balance equation, the temperature at customer j, given an unloading time t, can be estimated using Equation (3). Detailed derivation of Equation (3) is presented in Appendix A.

Temperature Prediction during Transportation
When the engine is running, the cooling equipment begins to remove the heat absorbed during the unloading process until the container temperature reaches Td. The rate of heat removal is denoted by Q c : m ij : Cargo mass between location i and j (kg) m a : Air mass (kg) s c : Specific heat of cargo (J kg −1 K −1 ) s a : Specific heat of the air (J kg −1 K −1 ) dT dt : Rate of change in temperature (K H −1 ) Applying the energy balance enables the prediction of container temperature during transportation between locations i and j (T ij ) using Equation (5). Detailed derivation of Equation (5) is presented in Appendix A. Equation (5) shows that the cooling equipment reduces the temperature at a rate of Q c (m ij s c + m a s a ) .
However, since the cooling equipment stops blowing cold air when the temperature reaches Td, the actual temperature is as follows:

CO 2 Emissions
The main source of CO 2 emissions in a perishable food distribution system is the fuel that is burned to provide energy for transport and refrigeration. The energy required for transportation, which is not specific to perishable products, depends primarily on the distance traveled, the weight of the vehicle and its cargo, the vehicle speed, and road/vehicle specifications. Bektaş & Laporte [39] developed a widely used method that integrates all these factors to predict road transportation energy consumption: w: Vehicle weight f ij : Weight of the load between node i to j v ij : Vehicle velocity d ij : Distance between nodes i and j α ij : Arc constant (i.e., road specification) β: Vehicle constant Bektaş & Laporte [39] calculated α ij and β as follows: a: Vehicle acceleration g: Gravitational constant C r : Rolling resistance θ ij : Slope of the road between locations i and j C d : Drag coefficient A: Vehicle frontal surface area ρ: Air density Stellingwerf et al. [17] added a method to calculate the energy consumed by refrigeration in temperature-controlled distribution. They illustrated that heat exchange between the air in the container and the ambient air during distribution is equal to the energy that the cooling system must remove to reduce the temperature inside the container. The MO-SVRP model presented in this paper extends their work by considering the air transfer ratio, air mass, unloading time, ambient air temperature, and predicted temperature inside the container (Section 3.1) to accurately estimate the heat exchange during unloading at location j using Equation (1). It is assumed that the heat exchange during transportation is negligible, since refrigerated containers are well-isolated. Considering both transportation and refrigeration energy consumption, the total energy consumption for a refrigerated vehicle v that is assigned to visit a set of locations L v is: x ij : Equals 1 if location j is visited immediately after location i; 0 otherwise Using the approach of Stellingwerf et al. [17], a refrigerated vehicle's total energy consumption can be converted to CO 2 emissions as follows: E v : Total CO 2 emissions of vehicle v (lb) TE v : Total energy consumption of vehicle v (kWh) µ: The efficiency of converting the chemical energy of the fuel to vehicle energy consumption (dimensionless) EC: Energy content of a gallon of fuel (kWh g −1 ) glb: Conversion factor for fuel: gallons to pounds (lb g −1 ) e: Conversion factor: fuel to emissions (dimensionless)

Food Product Freshness Based on Shelf Life Prediction
Increased temperature can elevate the growth rate of specific spoilage organisms (SSOs) in perishable food products [11]. Therefore, most shelf life prediction models require the temperature of the product over time, as well as product characteristics, to predict the remaining shelf life. Using the temperature prediction method explained in Section 3.1, an accurate estimate of product temperature from the time of vehicle departure from the depot until delivery is possible.
While most studies assume that the remaining shelf life of a product declines at a constant rate over time, the shelf life prediction model provided by Bruckner et al. [5] predicts a nonlinear increase in the number of SSOs in non-isothermal conditions. A product reaches the end of its shelf life when the number of SSOs reaches its maximum acceptable level. The product is not safe for consumption beyond this point and is considered spoiled. The remaining shelf life of food products can be estimated at any point in time, given the initial count of SSOs and characteristics of the food at presumed future storage temperature. Figure 2 shows the number of SSOs over time for a unit of product volume for constant temperature. : Conversion factor for fuel: gallons to pounds (lb g −1 ) : Conversion factor: fuel to emissions (dimensionless)

Food Product Freshness Based on Shelf Life Prediction
Increased temperature can elevate the growth rate of specific spoilage organisms (SSOs) in perishable food products [11]. Therefore, most shelf life prediction models require the temperature of the product over time, as well as product characteristics, to predict the remaining shelf life. Using the temperature prediction method explained in Section 3.1, an accurate estimate of product temperature from the time of vehicle departure from the depot until delivery is possible.
While most studies assume that the remaining shelf life of a product declines at a constant rate over time, the shelf life prediction model provided by Bruckner et al. [5] predicts a nonlinear increase in the number of SSOs in non-isothermal conditions. A product reaches the end of its shelf life when the number of SSOs reaches its maximum acceptable level. The product is not safe for consumption beyond this point and is considered spoiled. The remaining shelf life of food products can be estimated at any point in time, given the initial count of SSOs and characteristics of the food at presumed future storage temperature. Figure 2 shows the number of SSOs over time for a unit of product volume for constant temperature.  The Gompertz model [52] predicts the number of SSOs over time: N (t): SSO count (log10 cfu g −1 ) at time t, A: Initial SSO count of the food product at the time it is loaded into a refrigerated vehicle (log10 cfu g −1 ) C: The difference between the maximum SSO population level (a constant defined for each type of food product) and the initial SSO count A (log10 cfu g −1 ) M: Time at which the maximum growth rate is obtained (h) B: Relative growth rate at time M (h −1 ) The Gompertz model [52] predicts the number of SSOs over time: N (t): SSO count (log 10 cfu g −1 ) at time t, Sustainability 2020, 12, 6668

of 21
A: Initial SSO count of the food product at the time it is loaded into a refrigerated vehicle (log 10 cfu g −1 ) C: The difference between the maximum SSO population level (a constant defined for each type of food product) and the initial SSO count A (log 10 cfu g −1 ) M: Time at which the maximum growth rate is obtained (h) B: Relative growth rate at time M (h −1 ) The relative growth rate (B) is a function of temperature and is predicted by the Arrhenius equation [53]: F: Pre-exponential factor describing the number of times two molecules collide Ea: Activation energy for growth of SSOs (J mol −1 ) R: Gas constant (8.314 J mol −1 K −1 ) T: Absolute temperature (K) The freshness of the delivered products at location j (fr j ) can be evaluated using Equation (14) as the percentage of remaining shelf life before distribution to location j, given estimated time and temperature conditions. The remaining shelf life is the amount of time that it will take for the number of SSOs to climb from their current level to the maximum acceptable level.

Mathematical Model
The freshness prediction and CO 2 emission modules are integrated with a VRP model to create the MO-SVRP model. The perishable food distribution problem is defined as a directed graph in which a fleet of homogenous vehicles (V), each with a capacity of Q, deliver perishable food from a depot (node 0) to a set of customers (C) using a set of transportation paths (A) which connect the nodes. The order that fulfills the demand of customer i (d i ) must be delivered within the customer's required time window ([a i , b i ]). For ease of reference, all notations are given in Table 2.  The MO-SVRP model is formulated as a multi-objective mixed-integer program, and the mathematical formulation is as follows: Subject to The first objective (Equation (15)) minimizes transportation costs, including the cost of traveling between customer locations as the first term and dispatching costs associated with the fixed costs of using a vehicle in the distribution plan as the second term. The second objective (Equation (16)) maximizes the total freshness of the delivered products at each customer location. The third objective (Equation (17)) minimizes total CO 2 emissions generated by refrigerated vehicles during transit and unloading. Equation (18) ensures that each customer location can be visited by only one vehicle. Equation (19) prevents the load carried between two locations from being greater than the capacity of the vehicle. A vehicle can only leave the depot once (Equation (20)), and if a vehicle arrives at a location, it must also leave that location (Equation (21)). The amount of cargo unloaded at a customer location must equal that customer's demand (Equation (22)). Equation (23) ensures that a customer cannot be visited prior to the time that the previous customer is visited plus the unloading and travel times between these two customers. Similarly, Equation (24) prevents the first customer from being visited earlier than the time required to travel from the depot to that customer's location. Deliveries must occur within customers' required time windows (Equation (25)). A vehicle's load after leaving a customer location must be less than or equal to the capacity of the vehicle minus the demand of the visited customer (Equation (26)), and Equation (27) ensures that the load carried between two customer locations is at least equal to the demand of the next customer. Equations (28) and (29)

Solution Procedure
The MO-SVRP model presented in the previous section is difficult to solve. Even a VRP problem with a single objective and fewer parameters and variables is categorized as an NP-hard problem [24,54]. Consequently, a meta-heuristic approach was applied to solve the problem in a reasonable time. NSGA-II is an efficient and widely applied meta-heuristic approach introduced by Deb et al. [55] as a search technique for finding optimal solutions to multi-objective problems. This algorithm works based on iterative improvements in the pool of solutions' quality. In each iteration, genetic algorithm operators create offspring from the existing pool of solutions, and the solutions in this new pool are sorted based on the non-dominated sorting algorithm and crowding distance index. The top solutions are then selected for the next iteration. This section describes how this approach has been applied to solve the MO-SVRP.
In the NSGA-II algorithm, a "chromosome" represents a solution that assigns a list of customers to each vehicle in a particular delivery sequence. Each chromosome is an array consisting of n + v − 1 elements, in which n represents the number of customers and v represents the number of vehicles. An example is given in Figure 3. There are v − 1 special characters in the array ("*" in Figure 3), which divide the array into v sections (i.e., one section for each vehicle). The n remaining elements of the array are integer values from 1 to n, each of which is assigned to a customer. Thus, the sections between two special characters are lists of customers assigned to each vehicle, and the order of these numbers represents the delivery sequence. In the NSGA-II algorithm, a "chromosome" represents a solution that assigns a list of customers to each vehicle in a particular delivery sequence. Each chromosome is an array consisting of n + v − 1 elements, in which n represents the number of customers and v represents the number of vehicles. An example is given in Figure 3. There are v − 1 special characters in the array ("*" in Figure 3), which divide the array into v sections (i.e., one section for each vehicle). The n remaining elements of the array are integer values from 1 to n, each of which is assigned to a customer. Thus, the sections between two special characters are lists of customers assigned to each vehicle, and the order of these numbers represents the delivery sequence. 5   Initially, n customers and v − 1 special characters are randomly assigned to p chromosomes, where p is the size of the pool of solutions. Then, crossover and mutation operators are used to generate new solutions as the algorithm iterates.
The single point crossover operator is used to generate Pc chromosomes from the previous pool of solutions. Figure 4 provides an example, in which the first three elements the parent chromosomes are swapped to create two new offspring.  The crossover procedure is followed by a repair procedure [25], which is applied to fix chromosomes that have duplicated or missing customers or special characters. An example is shown in Figure 5.   Initially, n customers and v − 1 special characters are randomly assigned to p chromosomes, where p is the size of the pool of solutions. Then, crossover and mutation operators are used to generate new solutions as the algorithm iterates.
The single point crossover operator is used to generate P c chromosomes from the previous pool of solutions. Figure 4 provides an example, in which the first three elements the parent chromosomes are swapped to create two new offspring. In the NSGA-II algorithm, a "chromosome" represents a solution that assigns a list of customers to each vehicle in a particular delivery sequence. Each chromosome is an array consisting of n + v − 1 elements, in which n represents the number of customers and v represents the number of vehicles. An example is given in Figure 3. There are v − 1 special characters in the array ("*" in Figure 3), which divide the array into v sections (i.e., one section for each vehicle). The n remaining elements of the array are integer values from 1 to n, each of which is assigned to a customer. Thus, the sections between two special characters are lists of customers assigned to each vehicle, and the order of these numbers represents the delivery sequence. 5   Initially, n customers and v − 1 special characters are randomly assigned to p chromosomes, where p is the size of the pool of solutions. Then, crossover and mutation operators are used to generate new solutions as the algorithm iterates.
The single point crossover operator is used to generate Pc chromosomes from the previous pool of solutions. Figure 4 provides an example, in which the first three elements the parent chromosomes are swapped to create two new offspring.  The crossover procedure is followed by a repair procedure [25], which is applied to fix chromosomes that have duplicated or missing customers or special characters. An example is shown in Figure 5.  The crossover procedure is followed by a repair procedure [25], which is applied to fix chromosomes that have duplicated or missing customers or special characters. An example is shown in Figure 5. The crossover procedure is followed by a repair procedure [25], which is applied to fix chromosomes that have duplicated or missing customers or special characters. An example is shown in Figure 5.  Pm solutions are then randomly selected from the previous pool of solutions. The locations of two randomly selected elements from each chromosome are swapped, as shown in the example in Figure 6. Crowding distance is an estimate of the density of solutions around a particular solution in a front. The value of crowding distance for a particular solution is the summation of distances of the solutions with the neighboring solutions of the same front. Equations (31) and (32) show the mathematical equations by which crowding distance is calculated.
: Crowding distance of solution i D: Set of the problem domains In each iteration, binary tournament selection is applied to sort the solutions, first based on their ranks, and then based on their crowding distance.
The NSGA-II main loop consists of offspring generation and ranking and sorting modules.  Figure 7 shows how the pool of solutions in domains D1 and D2 is categorized in three fronts. Pm solutions are then randomly selected from the previous pool of solutions. The locations of two randomly selected elements from each chromosome are swapped, as shown in the example in Figure 6. Crowding distance is an estimate of the density of solutions around a particular solution in a front. The value of crowding distance for a particular solution is the summation of distances of the solutions with the neighboring solutions of the same front. Equations (31) and (32) show the mathematical equations by which crowding distance is calculated.
: Crowding distance of solution i D: Set of the problem domains In each iteration, binary tournament selection is applied to sort the solutions, first based on their ranks, and then based on their crowding distance.
The NSGA-II main loop consists of offspring generation and ranking and sorting modules. Crowding distance is an estimate of the density of solutions around a particular solution in a front. The value of crowding distance for a particular solution is the summation of distances of the solutions with the neighboring solutions of the same front. Equations (31) and (32) show the mathematical equations by which crowding distance is calculated. In each iteration, binary tournament selection is applied to sort the solutions, first based on their ranks, and then based on their crowding distance.
The NSGA-II main loop consists of offspring generation and ranking and sorting modules. Algorithm 1 presents a pseudocode that illustrates an iteration of the NSGA-II algorithm. At the end of each iteration, solutions with rank 1 are stored as Pareto-frontier. Algorithm 1 NSGA-II main loop pseudocode p = population of randomly generated chromosomes with size pops for i = 1 to Max number of iterations do for j =1 to (p c ÷ 2) do c 1 = 1 st randomly chosen parent chromosome c 2 = 2 nd randomly selected parent chromosome popc = append crossover (c 1 , c 2 ) for k = 1 to p m do m = a randomly chosen parent chromosome popm = append mutate (m) pop = merge (p, popc, popm) function (non-dominated sorting (input: pop)) return: Rank of chromosomes function (crowding distance (input: pop, rank)) return: crowding distance value for each chromosome function (sort population (input: pop, rank, crowding distance)) return: sorted pop based on 1) rank, 2) crowding distance pop = store only the top pop and truncate the others function (non-dominated sorting (input: pop)) return: Rank of chromosomes function (crowding distance (input: pop, rank)) return: crowding distance value for each chromosome function (sort population (input: pop, rank, crowding distance)) return: sorted pop based on 1) rank, 2) crowding distance Pareto_frontier = chromosomes with rank 1 Go to the next iteration if the stopping criteria are not met

Performance of the Solution Method
First, the performance of the NGSA-II solution algorithm in solving the MO-SVRP was tested on Solomon's datasets [56], which are widely applied to measure the quality of solutions for a VRP. The geographical distribution of the visiting location has a substantial impact on the performance of the VRP solution algorithm. So, these instances were provided in three categories: R, C, and RC, representing randomly generated, clustered, and mixed generated geographical data, respectively. Within these three categories, the MO-SVRP was solved for small instances (25 customers) and large instances (100 customers). However, the solutions should be compared with a competing algorithm to verify the efficiency and accuracy of the NGSA-II. Thus, weighted simulated annealing (w-SA) was used to solve the single objective weighted problem. According to the L1 metric method, the inverse of the optimum solution for each objective was used as the weight of the objectives in a single weighted objective function [23].
The NSGA-II and w-SA algorithms were coded in Python and run on a computer with 3.10 GHz Intel Core i9 CPU, 64 GB of RAM, and Windows 10 operating system. The parameter values that were used are given in Table 3. The algorithms terminate when the solutions in the Pareto front (for NSGA-II) or the optimum solution (for w-SA) do not improve after a certain number of iterations, which depends on the size of the problem. In Table 4, the column "w-SA" provides the values of each of the three objective functions for the best solution to each test problem instance. The values of the objective functions in the "NSGA-II" column correspond to the best solution that was found among the Pareto front solutions for each objective. The gap shows the difference between the objective values divided by the best value generated by these algorithms.
The results in Table 4 show that, on average, NSGA-II provides solutions that are 9.2%, 4.2%, and 8.0% better than w-SA in terms of cost, freshness, and emission objectives, respectively. The advantage of NSGA-II over w-SA is more pronounced for the cost objective when the size of the problem increases, or when the customers are more geographically clustered (i.e., the gap is largest for the C instances).

Optimality Analysis
Balancing cost, freshness, and emission objectives is necessary to improve the sustainability of perishable food distribution networks. The results in Table 5 demonstrate the conflicting nature of the three objectives: if the MO-SVRP problem is solved for a single objective, the values of the other two objective functions deviate significantly from optimality. Because the traveled distance is the primary driver of transportation cost and emissions, it is unsurprising that the emissions objective in the cost-optimal solution has only a 4% gap with the emission-optimal solution, and cost objective in the emission-optimal solution is only 14% higher than in the cost-optimal solution. However, keeping perishable products fresh requires faster delivery and fewer stops, which means that the vehicles may not be full. This increases the total traveled distance, which increases transportation cost and energy consumption, such that a solution that preserves 95% of the product's freshness results in 95% and 94% optimality gaps for cost and emissions objectives, respectively. In contrast, when a solution is cost-or emissions-optimal, freshness is 42% and 41% less than optimal, respectively. Most likely, there is more than one non-dominated solution in the Pareto frontier. So, a final solution that properly reflects the impact of all the objectives is chosen from the non-dominated set of solutions in the Pareto frontier. To choose a final solution, the approach developed by Bortolini et al. [32] was adapted. The presented method only ranks the solutions based on their distance to the optimum objective values, and other factors such as decision-maker priorities are not reflected in this method.
θ l : Represents a single calculated value for a solution l in Pareto frontier α l : Value of the first objective function for solution l α * l : Optimum value of the first objective function for solution l β l : Value of the second objective function for solution l β * l : Optimum value of the second objective function for solution l γ l : Value of the third objective function for solution l γ * l : Optimum value of the third objective function for solution l As shown in Table 5, the optimality gaps for each objective in the final solution are 52%, 57%, and 21% for cost, emissions, and freshness, respectively. When the only objective is to maximize freshness, the refrigerated vehicles are likely carrying loads that are far smaller than their capacity, rather than consolidating customer orders to ensure fast delivery and few delivery stops. In this scenario, the vehicles' traveled distances are very high, and consequently, distribution costs and CO 2 emission are at their highest level. Consolidating orders and increasing vehicle capacity utilization provides a substantial improvement in the cost and emission objectives, and the gaps of these objectives associated with the final solution reflect this.

Sensitivity to Shelf Life
To generate the results presented in the previous sections, it was assumed that the products have similar characteristics and the same shelf life (i.e., 2880 h). In reality, perishable products' shelf life can range from 168 h for highly perishable products, such as tomatoes, to 1440 h for moderately perishable products, such as oranges, and 2880 h for products with low perishability, such as apples. Therefore, the sensitivity of the MO-SVRP model to long, medium and short shelf life scenarios was analyzed.
The final values of the three objectives for these shelf life scenarios are shown in Table 6. The results emphasize the increase in cost and energy that is required to distribute more perishable items, as well as the loss of freshness for shorter shelf-life products. These results suggest that distributors should consider investing more money and time on improving container isolation and the efficiency of the diesel engine to reduce energy consumption and preserve quality when distributing products with a shorter shelf life. These results are compatible with the findings of Bortolini et al. [32] in which it was illustrated that products with lower shelf life have higher operating costs and produce a higher carbon footprint. Although there are recommended temperature ranges for perishable product storage, determining the specific temperature setting for a refrigerated shipping container can be challenging. Lower temperature settings preserve product quality, but this requires more energy. Therefore, the sensitivity of the MO-SVRP model to various temperature settings was assessed by solving the R101(25) instance at temperature settings between 263-269 degrees Kelvin (14-25 • F) in two-degree increments. Table 7 shows that the final solutions for the three sustainability objectives for different temperature settings inside the container. Increasing the temperature setting from 263 • K to 269 • K causes a 15% reduction in product freshness (from 75% to 60%). Because the MO-SVRP recommends faster delivery to compensate for the loss in quality that results from an increase in temperature, the transportation distance and cost tend to increase at higher temperature settings. However, the effect of an increased temperature setting on emissions is more complex. On one hand, transportation consumes more energy due to the increase in traveled distance. On the other hand, energy consumption for refrigeration decreases. This leads to a decrease in overall energy consumption throughout distribution.

Conclusions and Future Research
In this paper, the VRP is extended to consider the perishability of the food products and refrigerated vehicles' CO 2 emissions in a multi-objective framework. The model addresses sustainability concerns associated with the perishable food distribution problem in a mathematical model by minimizing transportation costs and CO 2 emissions and maximizing product freshness. Integrating a temperature estimation model provides accurate predictions of refrigeration-related emissions and product quality loss during distribution. The presented MO-SVRP model was solved using an adapted NSGA-II algorithm, and the performance of the solution algorithm was tested and verified against the w-SA algorithm over the Solomon [56] data sets.
The results of the analysis revealed that the sustainability objectives are conflicting, such that optimizing one objective can degrade the other two objectives. In particular, optimizing the freshness objective leads to a solution in which deploying all the available vehicles to minimize delivery time and the quality-deteriorating effect of temperature abuses at the delivery locations. Therefore, the optimality gaps in distribution costs and CO 2 emission objectives are high when only the freshness objective is optimized. Furthermore, when vehicles are carrying highly perishable products, the MO-SVRP model recommends faster delivery and fewer stops, which increases energy consumption and distribution costs. Sensitivity analysis over different temperature settings inside the shipping container indicates that even small increments in the temperature settings can have a huge impact on CO 2 emissions and freshness objectives. This sensitivity analysis can be a helpful tool to determine the best temperature setting to achieve sustainability goals.
This paper is mainly focused on presenting a methodology to accurately measure and optimize sustainability goals in the distribution of perishable food products. It is highly recommended to apply the methodology in practice to evaluate the sustainability goals with respect to the real-world parameters and constraints. The methodology of this research is limited to the single product distribution system with single compartment refrigerated vehicles. Considering the growth in the application of multi-compartment refrigerated vehicles to distribute multiple types of products with one vehicle, it would an interesting expansion for this study to measure and optimize sustainability objectives for perishable food distribution with multi-compartment refrigerated vehicles in which each food product has different perishability parameters and the compartments of the refrigerated vehicles can be set to different temperatures.
The outcomes of this research highlight the necessity of integrating and accurately estimating multiple influential factors that impact sustainability in a perishable food distribution network to find solutions that are cost-effective, reduce food waste, and decrease emissions generated by refrigerated vehicles.