THERMAL PERFORMANCE OF REFRIGERATED VEHICLES IN THE DISTRIBUTION OF PERISHABLE FOOD

Abstract

The temperature of refrigerated products along the distribution process must be kept within close limits to ensure optimum food safety levels and high product quality. The variation of product temperature along the vehicle routing sequence is represented by non-linear functions. The temperature variability is also correlated with the time required for the refrigerated unit to recover after cargo unloading, due to the cargo discharging process. The vehicle routing optimization methods employed in traditional cargo distribution problems are generally based on the Travelling Salesman Problem with the objective of minimizing travelled distance or time. The thermal quality of routing alternatives is evaluated in this analysis with Process Capability Indices (PCI). Since temperature does not vary linearly with time, a Simulated Annealing algorithm was developed to get the optimal solution in which the minimum vehicle traveling distance is searched, but respecting the quality level expressed by a required minimum PCI value.

cold chain; vehicle routing; perishable cargo distribution; PCI; TTI

1 INTRODUCTION

Lifestyle changes over the past decades led to increasing consumption of refrigerated and frozen foods, which are easier and quicker to prepare than the traditional types of food. In order to ensure product quality and health safety (^{Andrée et al., 2010}; ^{Wegener, 2010}; ^{Daelman et al., 2013}), the control of temperature throughout the cold chain is necessary. In fact, a number of factors affect the maintenance of quality and the incidence of losses in fresh food products, such as (a) the initial quality of the commodity; (b) the temperature at which the product is held during handling, storage, transport, and distribution; (c) the relative humidity of the postproduction environment; (d) the use of controlled or modified atmospheres during storage or transit; (e) chemical treatments for the control of decay or physiological disorders; (j) heat treatments for decay control; and (g) packaging and handling systems (^{Harvey, 1978}). But, since temperature largely determines the rate of microbial activity, which is the main cause of spoilage of most fresh food products (^{Tarantilis & Kiranoudis, 2002}; ^{Mendoza et al., 2004}; ^{Giannakourou et al., 2005}; ^{James & James, 2010}; ^{Ellouze &Augustin, 2010}; ^{Mai et al., 2011}), continuous monitoring of the full time temperature history usually allows for an adequate control of the process along the short and medium distance distribution situations (^{Giannakourou et al., 2005}). The quality of these products might change rapidly because they are submitted to a variety of risks during transport and storage that are responsible for material quality losses (^{Pereira et al., 2010}; ^{Tassou et al., 2012}). Metabolic activities generally increase as storing temperatures are elevated (^{Borch & Arinder, 2002}; ^{McMeekin et al., 2008}; ^{Havelaar et al., 2010}; ^{Mai et al., 2011}). On the other hand, short interruptions in the control of the cold chain may lead to quick deterioration of product quality (^{Garcia, 2008}).

Therefore, the required product temperature range should be maintained from production to consumption. However, the main difficulties are encountered at the weakest links or interfaces of the cold chain as delivering, loading or unloading operations and temporary storage where products are generally handled in somewhat uncontrolled temperature ambiences. For an adequate control of the cold chain, some practical recommendations have been presented in the literature. Many investigations have been carried out over the past years on the logistics food chain in order to characterize the temperature variation as a function of product properties, ambient conditions, the size and kind of packaging, the presence and thickness of air layer between the product and the package, solar exposure, and the use of insulating pallet covers. These studies show large variations of temperature in practice (^{Moureh & Derens, 2000}).

In practical terms, the maintenance of an optimal temperature throughout the postproduction handling chain is one of the most difficult tasks, and is far from being universally attained. Even when road transport can provide satisfactory temperatures within the limits of acceptability, the displacement time may be too long for short-life products to be transported over long distances. On the other hand, in short or medium distance delivery runs, the chilled or frozen product can be subjected to many door openings, where there is heat ingress directly from outside air and from personnel entering to select and remove product (^{James et al., 2006}; ^{Pereira et al., 2010}). Frequent door openings can also lead to increased evaporator frosting, resulting in a reduction of the evaporator's performance and an increase in the need for defrosts, particularly in humid weather conditions (^{Estrada-Flores & Eddy, 2006}). Additionally, the design of the vehicle refrigeration system has to allow for extensive variation in load distribution, which is a function of different delivery rounds, days of the week and the removal of product during a delivery run. A refrigeration system's ability to respond to sudden demands for increased refrigeration is often restricted by the power available from the vehicle. As a result, there are substantial difficulties in maintaining the temperature of refrigerated products transported in small and medium-size refrigerated vehicles that perform multi-drop deliveries to retail stores and caterers (^{James et al., 2006}).

This paper reports a Time-Temperature Indicators (TTI) analysis of a distribution of refrigerated products along a route containing a number of retail customers with different demand levels. Time-Temperature Indicators (TTI) can effectively be used to monitor the time/temperature impacts on product quality, offering a cost-effective way to detect problematic points in the chill chain. TTI are defined as specific elements collected by electronic equipment that show easily measurable, time and temperature dependent changes that cumulatively indicate the thermal history of the product from the point of manufacture to its destination (^{Giannakourou et al., 2005}; ^{Estrada-Flores & Eddy, 2006}; ^{Sahin et al., 2007}). One simple type of TTI application involves the registration of temperature measurements along a pre-defined time horizon and the subsequent analysis of this historical data to get the corresponding probabilistic temperature distribution, its average and standard deviation, unexpected temperature rises and drops, and other elements which occur along the transport, storage, and distribution processes (^{Simpson et al., 2012}; ^{Pereira et al., 2010}; ^{Estrada-Flores & Eddy, 2006}). Process Capability Indices (PCI), on the other hand, can additionally be calculated to yield easily computed evaluation coefficients measured with dimensionless functions on TTI parameters and specifications (^{Chang & Bai, 2001}; ^{Chang et al., 2002}; ^{Estrada-Flores & Eddy, 2006}). This kind of TTI application helps to reveal undesired thermal conditions that may impair the compliance of product quality requirements along the cold chain.

The paper analyses alternative vehicle routing strategies intended to minimize travel cost, but at same time keeping thermal PCI performance indicators within the required levels. It is shown that the standard TSP (Travelling Salesman Problem) approach, used to solve classical routing problems where vehicle travel distance or time isminimized, usually leads to temperature restriction violations. Thus, instead of using a heuristic to get the optimized TSP vehicle routing sequence, such as the largely employed 2-opt and 3-opt improvement methods (Syslo et al., 2006), other routing strategies, apart from the TSP criterion, can be employed (^{Ferreira & Pureza, 2012}; ^{Faraz et al., 2013}). In this paper a Simulated Annealing algorithm was specifically developed to optimize the routing problem.

2 THERMAL ASPECTS OF REFRIGERATED TRANSPORT

Thermal PCI serve as benchmarking elements to compare refrigerated transport systems and to ensure that these systems provide a minimum level of operating effectiveness (^{Estrada-Flores & Eddy, 2006}). Refrigerated freight transport must comply with international ISO rules (^{ISO, 1996}). In addition, Europe, the United States and many countries have their own specific rules. In Brazil, there is a specific regulation on the subject issued by ABNT - Brazilian Association of Technical Norms (^{NBR 14701, 2001}). EU legislation covers temperature control requirements during the storage and transport of perishable foods. In Europe, the European Commission revised the related rules in 2006 and regulation EC Nº 852/2004, on the Hygiene of Foodstuffs, requires manufacturers to have suitable temperature controlled handling and storage facilities that are able to maintain food at appropriate temperatures, also permitting that these temperatures be adequately monitored, controlled and recorded (^{Tassou et al., 2012}). The Australian regulation requires that refrigerated vehicles have a sufficient cooling capacity to offset the heat load infiltrated through the insulated body, plus an extra capacity calculated as a percentage of the heat leakage. For trucks submitted to only one door opening during an eight-hour period, the percentage extra cooling capacity required is 25%. The extra capacity increases to 200%, for trucks that have 31 to 35 door openings per day. The reserve power is also related to the time-based deterioration of both insulation and refrigerating plant, the differences between testing and in-service conditions, the heat of respiration while transporting fresh produce and the extra-cooling capacity necessary to eliminate defrosting heat (^{Estrada-Flores & Tanner, 2005}). In local delivery, which is the focus of this paper, there is also intensifying demand from legislation and retailers for lower delivery temperatures, and increased pressure on fleet operators to improve temperature control (^{James et al., 2006}).

Road transport refrigeration equipment is required to operate reliably in much harsher environments than stationary refrigeration equipment. Due to the wide range of operating conditions and constraints imposed by available space and weight, transport refrigeration equipment has lower efficiency levels than stationary systems. The thermal performance of refrigerated vehicles is also dependent of age, since insulation materials deteriorates with time due to the inherent foam characteristics. Recent data show a typical loss of insulation value of between 3% and 5% per year which can lead to considerable rise in the thermal conductivity after a few years. If a 5% yearly ageing is assumed, a vehicle after nine years of operation may show an increase of approximately 50% in energy consumption and CO₂ emissions (^{Tassou et al., 2009}, ²⁰¹²).

The effects of door openings are another important concern in refrigerated cargo transport (^{Estrada-Flores & Eddy, 2006}; ^{Pereira et al., 2010}). ^{Tso et al. (2002)} used a commercial computational fluid dynamics (CFD) program to model the effects of door openings on air temperature within a refrigerated truck. They carried out a series of experiments to study the effect of door openings with unprotected doors, with air curtains and with plastic strip curtains, showing expressive temperature raises during the unprotected operation.

3 THERMAL PERFORMANCE ANALYSIS

Generally, models that address the prediction of heat and mass transfer during transport can be divided into those that consider the environment within the transport unit (usually in regard to the airflow) and those that concentrate on the temperature of the product (^{James et al., 2006}). Some models combine these aspects and deal with the temporal aspects of transportation: fluctuating ambient conditions, door openings, product removal/loading, etc. Other models specifically address the effects of transportation temperatures on microbial growth and its influence on food safety.

With regard to models based on the ambience in refrigerated vehicles, ^{James et al. (2006)} resumed a number of approaches that started in the 1990s. These models use computational fluid dynamics (CFD) techniques to analyse airflow in cold environments and through cold-store doorways (^{Moureh & Derens, 2000}; ^{Ge & Tassou, 2001}; ^{Moureh et al., 2002}). The effects of door openings are an important concern in refrigerated cargo transport (^{Estrada-Flores & Eddy, 2006}; ^{Pereira et al., 2010}). ^{Moureh & Derens (2000)} employed a CFD model to study temperature rises in pallet loads of frozen food during distribution. They specifically looked at the times during loading, unloading and temporary storage when the pallets would be in an environment above 0ºC.

A number of papers on refrigerated food transport address the time-temperature evaluation along the cold chain process (^{Gigiel et al., 1998}; ^{James & Scholfield, 1998}; ^{Tsironi et al., 2001}; ^{Smolander et al., 2004}; ^{Giannakourou et al., 2005}; ^{Estrada-Flores & Eddy, 2006}; ^{Pereira et al., 2010}; ^{Nga, 2010}). The main purpose of maintaining good temperature control during refrigerated transport is to decrease the rate of microbial growth and hence maintaining the safety and eating quality of the food (^{Zhou et al., 2010}). In fact, there are many microbial models that can be applied to modelling the growth of microorganisms in food during transport. But, according to ^{James et al. (2006)}, relatively few studies have been published specifically on the refrigerated transport subject. Even fewer of these have looked at the fully integrated approach of combining dynamic microbial growth modelling with heat and mass transfer models that could represent the characteristics of a food distribution system (^{James et al., 2006}). However, in shortdistance transportation of refrigerated food products, as in the case of this work, it suffices to observe that the temperature of the product is maintained within pre-established limits. This condition is adopted in a number of studies reported in the literature based on Time-Temperature Indicators (TTI).

Papers investigating the thermal performance of refrigerated vehicles can be classified into four groups: (a) approaches of pure theoretical nature, (b) performance laboratory tests (c) field data gathering, and (d) TTI based simulations. Studies of the first group use pure mathematical CFD models based on physical aspects of heat transfer (^{Cuesta et al., 1990}; ^{Zhang et al., 1994}, ^{Campañone et al., 2002}; ^{Flick et al., 2012}; ^{Hoang et al., 2012a},^b). In the second group are the papers describing controlled laboratory tests, such as ^{Moureh & Derens (2000)}, ^{Tso et al. (2002)}, ^{Estrada-Flores & Eddy (2006)}, and ^{Garcia (2008)}. In the third group are the papers involving field tests, with information eventually supplemented with laboratory data, such as the CoolVan experiment (^{Gigiel et al., 1998}; ^{James & Scholfield, 1998}; ^{James et al., 2006}), as well as other similar efforts (^{Giannakourou et al., 2005}; ^{Pereira et al., 2010}; ^{Nga, 2010}). The fourth group comprises TTI - Time-Temperature Indicator analyses, usually employing computer simulations (^{CoolVan Manual, 2000}, ^{Hoang et al., 2012a},^b). Of course, a good part of these papers employ combined methods. In this paper, the approach (d) will be used, involving time-temperature indicators (TTI) associated with simulation of the thermal performance of refrigerated vehicle tours.

For the purposes of this study, that focusses a practical logistics application, the employment of pure theoretical models are not recommended due to the great number of variables and situations, and the dynamic nature of the problem. The performance of thermal laboratory tests, by its turn, would allow for the gathering of more reliable technical data, but would not include important in-the-field behavioural characteristics and information. Besides, the analytical methodology presented in this work is essentially intended to be used in the planning phase of the problem, which requires an approximate accuracy level. For further applications to real problems, more accurate data should be sought, and laboratory tests could eventually be contemplated in such situations. Field tests to gather thermal data, on the other hand, would be acceptable in developing countries like Brazil only under rigid technological and experimental control, a condition hard to achieve presently (^{Pereira et al., 2010}). As a result, approach (d) was selected for this application, involving time-temperature indicators (TTI) associated with thermal performance simulation of refrigerated vehicle tours.

A number of papers on refrigerated food address the TTI evaluation along the cold chain process (^{Gigiel et al., 1998}; ^{James & Scholfield, 1998}; ^{Jacxsens et al., 2002}; ^{Giannakourou & Taoukis, 2003}; ^{Giannakourou et al., 2005}; ^{Estrada-Flores & Eddy, 2006}). The main purpose of maintaining good temperature control during refrigerated transport is to decrease the rate of microbial growth and hence maintaining the safety and eating quality of the product. In fact, there are many microbial models that can be applied to represent the growth of microorganisms in food during transport. But since temperature largely determines the rate of microbial activity, which is the main cause of spoilage of most fresh and frozen food products, continuous monitoring of the full time temperature history usually allows for the adequate control of the process along short and medium distance distribution situations.

One of the most systematic attempts to predict the temperature of refrigerated food during multidrop deliveries has been the CoolVan research programme. CoolVan is a software developed by the Food Refrigeration and Process Engineering Research Centre at the University of Bristol, UK. The brief description set forth was extracted from ^{Gigiel (1997)}, ^{Gigiel et al (1998)}, and from the ^{Cool Van Manual (2000)}. The objective of CoolVan is to aid the design and operation of small and medium delivery vehicles intended to distribute refrigerated food products (^{Gigiel, 1997}; ^{Gigiel et al., 1998}; ^{James & Schofield, 1998}; ^{James et al., 2006}; ^{CoolVan Manual, 2000}). The software contains a mathematical model that predicts food temperatures inside a refrigerated delivery vehicle, analysing the temperature changes that take place during a delivery journey as well as the energy used by the refrigeration equipment. The model is solved using an implicit finite difference method. It starts with the given initial conditions and proceeds to the end of the journey with variable time steps. The time step is increased when there are few events during the journey. However, the time step is reduced and the program slows down in order to model in detail any event which occurs along the route as, for example, when the vehicle stops and the door is opened. In this way the program is able to model the details of the journey accurately, while at the same time, running at the fastest possible speed. It shows the thermal consequences of changes in delivery patterns, theweather, the vehicle construction characteristics, its refrigeration system, and the type and mass of the food transported.

There are substantial difficulties in maintaining the temperature of chilled foods transported in refrigerated vehicles that conduct multi-drop deliveries to retail stores and caterers. The vehicles have to carry a wide range of products and operate under diverse ambient conditions. During any delivery run, the chilled product can be subjected to many door openings, where there is heat ingress directly from outside and from personnel entering to select and remove product. On the other hand, the design of the refrigeration system has to allow for extensive differences in load distribution, dependent on different delivery rounds, days of the week and the removal of product during a delivery run. A refrigeration system's ability to respond to sudden demands for increased refrigeration is often restricted by the power available from the vehicle. All these problems converge to produce a complex interactive system. With so many interacting variables it is not feasible to obtain, by purely theoretical means, the reasoning required to design and to operate vehicles that wouldmaintain the required food temperatures. Thus, the CoolVan research project was set up to provide a predictive model that would assist in specifying the design of, and the equipment for small and medium delivery vehicles. Operators of the vans need to know in advance whether a particular vehicle, on a particular route, under given ambient conditions, will be able to deliver food at the correct temperature levels (^{Gigiel et al., 1998}).

The heart of the CoolVan program is the temperature of air inside the vehicle. The internal air exchanges heat with the outside environment by the movement of air into and out of the truck, while the doors are either opened or closed. Heat is conducted into the van through the insulation from the ambient air and from solar radiation on the outside of the van. Inside the van, the racking, fittings, trays and food exchange heat with the van air. Evaporator coils or eutectic plates cool the inside air to maintain its temperature. The composition of each side of the van is modelled separately and heat flow through each side assessed independently. Each side can be modelled as up to three layers, an inner and outer skin and insulation. The program models the heat transfer through each wall of the vehicle separately. However, the walls of vehicles are not uniform so the heat flowing through the walls will vary. Average values are then calculated for the thermal properties. An important thermal effect is due to air infiltration into the vehicle through gaps and cracks as it moves along the road network. It is very difficult to know how much air enters into the vehicle and, indeed, it will vary after each instance that the door has been opened and closed. These detailed measurements cannot be accurately modelled with CoolVan software, but its results are accurate enough for planning purposes and approximate operation analyses.

The usual food distribution scheme starts with the vehicle being loaded at the distributor's premises and travelling next to a series of retail outlets, where the individual lots are discharged in sequence. Often, the vehicle has a large number of servicing stops in a journey, when the doors are opened and food is removed. Sometimes, food which has passed its shelf-life date, together with empty trays, return from the retail shops to the distributor. In another cases the insulation, door protection and refrigeration equipment fitted to the vans prove to be inadequate to maintain food temperature as cold as required. The food is assumed to be at the same temperature anywhere in the vehicle along the journey. Additionally, it is assumed that the food is loaded into the vehicle at the same uniform temperature. In practice, this is rarely the case. There are temperature gradients within each item of food and different foods will have different temperatures. But, as a first approximation for logistics analysis, such as the one in this paper, the even temperature approximation is satisfactory.

The CoolVan program is divided into three parts. The first part allows the user to change over a large number of different variables that are input into the main program: the weather, the way the vehicle is used, its construction, the refrigeration plant and the product that is put in it. The second part of the program does the actual calculations and puts these into an output file that can be saved for future reference or viewed in another application. The last part of the program plots the output data graphically on the screen to allow the user to see what is going on during the journey. Vehicle data are fed into the CoolVan program: the thermal properties of the insulation system, the year of the van manufacture, the ageing rate which depends on the vehicle maintenance characteristics, etc. Then, the program calculates the reduced thermal properties of the vehicle insulation. The mathematical structure of the program also allows different external heat transfer coefficients to be entered for each side of the vehicle. Solar radiation onto each surface of the van is modelled separately. The infiltration of outside air into the van is dependent on the van structure, the degree of maintenance and the speed of the vehicle. These effects were measured empirically in several vans, allowing for the fitting of appropriate equations and parameters into the model.

Controlling temperature during Infiltration of air into the vehicle during door openings is also an important issue. Once the door is opened, the amount of air infiltrating the vehicle approximately increases linearly with time (^{Gigiel et al., 1998}). Vertical air curtains, plastic strip curtains, inner lightweight sliding doors or shielding are sometimes used on doors in order to limit the air exchange through them. To accurately model the air exchange in a particular van it is necessary to perform empirical measurements. However, real measurements of the air exchange rate on several vans showed that a reasonable approximation is that it varies as a function of the height and width of the door only.

4 PROCESS CAPABILITY INDICES TO ASSESS TTI DATA

A process capability index (PCI) is a numerical element that compares the characteristics of a production or servicing process to engineering specifications (^{Pearn & Chen, 1999}; ^{Chang & Bai, 2001}; ^{Chang et al., 2002}; ^{Czarski, 2008}). A value of such an index equal or larger than a pre-established level indicates that the current process is capable of producing results that, in all likelihood, will meet or exceed the pre-defined requirements. Process Capability Indices (PCI) are frequently used as an integral part of the statistical control of process quality and productivity. A capability index of this sort is convenient because it reduces complex information about the quality of the process to a single number. Thus, process capability indices can be defined as measurements of the ability of the production and servicing processes to meet the specifications. In our application, the analysis of the temperature variability inside a refrigerated vehicle is performed by means of PCI estimation applied to TTI values, showing the temperature distribution on the data collected along a typical distribution journey. We closely follow ^{Chang & Bai (2001)}, ^{Chang et al. (2002)}, ^{Czarski (2008)} and ^{Estrada-Flores & Eddy (2006)} in applying PCI analysis to the distribution of refrigerated food products.

4.1 Capability Indices for Normally Distributed Data

Usually, capability indices are employed to relate the process parameters to engineering specifications that may include unilateral or bilateral tolerances, with or without a target value (nominal value). The resulting indices are dimensionless and provide a common, easily understood way for quantifying the performance of a process. In this application, the monitored variable is the temperature θ inside the vehicle along a typical distributionjourney of refrigerated food products, with mean μ and standard deviation σ. In this case there are two-sided specification limits for θ, respectively the upper value USL and the lower value LSL. Four capability indices are commonly used for variables normally distributed: C _p, C _pk, C _pm and C _pmk (^{Gonçalez & Werner, 2009}). The C _p index is defined as

Clearly, the aim of process control is to make C _p as large as possible. The adoption of the C _p index presupposes that the variable θ is normally distributed and that the mean μ is equal to the specified target value θ _T. The coefficient C _p is defined as the ratio of the allowable tolerance spread and the actual spread of the data. If C _p > 1, it indicates that the temperature variation fits within the specified limits. In practice, if C _p < 1, the process does not meet the defined specifications; if 1 ≤ C _p < 1.33 the process probably meets the requirements but additional attention must be taken; if C _p ≥ 1.33 the process is fully capable (^{Gonçalez & Werner, 2009}). The six-sigma coverage represents the spread of 99.73% of the data in normally distributed processes. On the other hand, the eight-sigma coverage, which represents the case with C _p ≥ 1.33, covers practically 100% of the data in similar circunstances. In practice, there are cases in which the process is not centred on the target value θ _T , i.e μ ≠ θ _T . To avoid such drawback, the index C _pk is defined

On the other hand, the index C _pm explicitly considers the distance between the mean and the target value θ _T (^{Gonçalez &Werner, 2009}):

For practical purposes, sample estimators can be used for calculating the above indices.

4.2 Capability Indices for Non-Normal Data

^{Gonçalez & Werner (2009)} compared a number of approximate methods to evaluate non-normal situations and suggested that the method of ^{Chen & Ding (2001)} reflects with accuracy the number of non-conforming items in the evaluated sample, being superior to other methods. Such method adjusts the values of PCIs in accordance to the degree of skewness of the underlying population, by using specific factors in computing the deviations above and below the variable mean. The method is based on the idea that σ can be divided into upper and lower deviations, σ _U and σ _L , which represent the dispersions of the upper and lower sides around the mean μ, respectively. An asymmetric probability density function f (θ) can be approximated with two normal pdfs (^{Chen & Ding, 2001})

with the same mean μ but different standard deviations 2σ _U and 2σ _L, where ∅ represents the standard normal pdf. The upper and lower sides of f (θ) are approximated with the upper side of f _U (θ) and the lower side of f _L (θ), respectively. The values of σ _U and σ _L are computed as (^{Chang et al., 2002})

The value of C _p is defined as (^{Chang & Bai, 2001})

Making D _θ = 1 + |1 − 2P _θ |, expression (1) can be simplified to

where 1/D _θ is a corrective coefficient on (1) due to the skewness of the probability distribution of θ. On the other hand, according to ^{Chang et al. (2002)}, the value of C _pk corrected for skewness can be estimated as follows. First, the upper and lower capability indices are defined as

The PCI of the temperature θ, for the case of two-side specifications (upper and lower values) is

In Equations (8) and (9), 2σ _U and 2σ _L are used in place of σ to reflect the degree of skewness in the probability distribution of θ. When the underlying distribution of θ is symmetric, P _θ = 0.5, repeating the value given by (1). However, if it is skewed, the value of C _pk given by (10) is smaller than the standard one given by (1). The numerical estimation of C _p and C _pk is performed as follows. Let (θ ₁ , θ ₂ , . . . , θ _n ) be the sample containing n values of the temperature θ. The mean and the standard deviation of the sample are computed. The probability P _θ can be estimated by using the number of observations less than or equal to (^{Chang et al., 2002}):

where g(x) = 1 for x ≥ 0 and g(x) = 0 for x < 0. Then C _pk can be estimated by substituting , , and _θ for μ, σ and P _θ respectively in Eq. (10). In this application the series of TTI data are skewed (see Section 6.4), and the method of ^{Chen & Ding (2001)} has been employed accordingly to get C _pk values.

A fourth capability index C _pmk is used when the mean μ departs from the target value θ _T (^{Gonçalez & Werner, 2009}). In this application there is no target value for the temperature θ, being only necessary to respect the two-sided specified limits, namely LSL and USL. Consequently, the index C _pk was adopted.

5 PROBLEM DESCRIPTION

The objective of the analysis is to investigate a regional distribution of ready-to-eat refrigerated meat products (ham, turkey and chicken breasts, salami, sausage). Other studies analysing distribution schemes of perishable products are ^{Oswald & Stirn (2008)}, ^{Yu & Nagurney (2013)}, and ^{Faraz et al. (2013)}. The urban distribution district is located about 84 km from the base depot. The served urban district has an approximated area of 73 sq.km, where 12 retail shops are located, as shown in Figure 1. Traditionally, the optimal sequence of points to be visited is obtained with a TSP (Travelling Salesman Problem) algorithm, that yields the shortest Hamiltonian cycle, which includes, in this application, all the retail outlets, plus the depot. Figure 1 depicts the TSP route obtained with a 3-opt local search heuristic (^{Syslo et al., 2006}), with a total extension of 204.1 km.

Figure 1
Distribution of refrigerated food in the district: the TSP vehicle route.

Along the delivery journey, the air temperature inside the truck tends to increase when it stops to discharge the products, and the temperature inversely decreases while the vehicle is travelling from a client to the next one in the routing sequence. Normally, a heavier discharge takes more time to be accomplished. If the next visiting point is located close by, the corresponding travelling time might be too short, and the temperature would not decrease enough to compensate for the heat surge observed in the last unloading stops. This means that the vehicle routing sequence has an important influence on the TTI history, indicating that both have to be investigated jointly. When the cargo mass inside the vehicle decreases, the temperature tends to decrease sharply during the local travel segments. In the CoolVan simulator, a thermostat turns off the refrigeration system when the temperature reaches 2◦C, in order to avoid chilling injuries. Conversely, during local travel segments, the refrigeration equipment turns on when θ 5◦C, but this does not happen during cargo unloading, since the refrigeration module and the truck engine are not in operation (to avoid discharging gases to enter the cargo compartment). This situation is shown in Figure 2 when the truck is travelling from client 10 to client 11, and from client 12 to client 3.

Figure 2
Graphical representation of TTI results obtained with the CoolVan software for the TSP delivery route.

Normally, the search for an optimal vehicle routing sequence requires quite a number of combinatory evaluations. On the other hand, the process of obtaining accurate TTI data is not a simple task since it requires special laboratory settings and/or elaborate field tests. One alternative possibility is to apply the simulation approach, such as the CoolVan software (Section 3), in order to gather basic TTI data to be used in a computer aided routing analysis. This is because the CoolVan program does not permit automatic replications of the runs, whereas combinatorial routeing analysis implies automatic combinatorial changes of the delivering sequences in order to obtain the resulting C _pk coefficients and the traveling distances.

The model employed to solve the proposed problem involves a combination of methods divided in four steps. First, the CoolVan software has been used to simulate the thermal characteristics of the basic routing schemes, one of them represented by the TSP formulation as shown in Figure 1. Next, taking the CoolVan simulation results, the operating stages that compose the temperature evolution along the distribution journey are analysed individuallyin order to define mathematical functions that relate temperature to the explaining variables. Third, a computer program has been developed to estimate TTI values, step by step, along each routeing sequence, employing the mathematical functions adjusted in the last step. A PCI coefficient is computed for each route scheme in order to analyse the thermal performance of the alternative delivery sequences. Next, and associated with this, a Simulated Annealing algorithm (Section 6) analyses each candidate route, searching for the minimum travel distance configuration, but at same time respecting the restriction C _pk ≥ 1.33. The methodology will be described in more detail in the next sections.

This line of research has some points in common with the paper by ^{Oswald & Stirn (2008)}. Their approach differs from the present investigation in two aspects. First, the variable that expresses product quality in our research is temperature, whereas in ^{Oswald & Stirn paper (2008)} product quality (fresh vegetables) is based on market acceptance: one has 100% qualitywhen the product can be sold entirely at the current market price and the quality drops to 0% when the product loses completely its commercial value (^{Oswald & Stirn, 2008}). Second, themathematical routing model developed by these authors builds up the routes to be assigned to the vehicles step by step, while in this application there is only one truck and the route sequence is defined in a combinatorial way.

5.1 CoolVan route simulation

The CoolVan software is still not available commercially but its developers kindly accepted to ran basic routing configurations to serve as a data basis for this application. The necessary inputs for the CoolVan simulation are presented in Table 1. Next, the data of the food products are fed into the CoolVan program (^{CoolVan Manual, 2000}). The product is ready-to-eat refrigerated meat products (ham, turkey and chicken breasts, salami, sausage). A total of 12,000 kg of assorted products are distributed in the daily round, with two retailers receiving larger quantities (7,000 and 2,000 kg respectively), while the other ten clients getting 300 kg each (Table 2). The basic delivery information shown in Table 2, and in Figure 2, illustrates the TSP route. Table 2 also shows the distances travelled along the various segments of the TSP route and the delivering times at the retail premises. It is worth to notice, in Figure 2, the sharp changes of the internal vehicle temperature in the last phases of the simulation process, when the product mass inside the van gets successively smaller.

Table 1
CoolVan simulation inputs.

Table 2
Input data for the TSP sequence of delivery points.

Stage type	Origin	Destination	Travel distance (km)	Travel time (min)	Product delivered (kg)	Unloading time (min)
Line-haul	Depot	Client #1	84.3	85.6	300	10
Local travel	Client #1	Client #4	1.6	3.1	300	10
Local travel	Client #4	Client #5	2.9	5.8	300	10
Local travel	Client #5	Client #8	1.2	2.3	7,000	64
Local travel	Client #8	Client #7	0.9	1.9	2,000	24
Local travel	Client #7	Client #6	1.3	2.6	300	10
Local travel	Client #6	Client #9	2.5	5.1	300	10
Local travel	Client #9	Client #10	4.2	8.3	300	10
Local travel	Client #10	Client #11	5.7	11.3	300	10
Local travel	Client #11	Client #12	0.6	1.2	300	10
Local travel	Client #12	Client #3	6.2	12.4	300	10
Local travel	Client #3	Client #2	3.8	7.5	300	10
Line-haul	Client #2	Depot	86.4	85.6	-	-
Total			204.1	232.7	12,000	188

The results of the CoolVan simulations are saved in Excel format. Figure 2 exhibits the TTI results obtained with the CoolVan software for the TSP sequence. The simulation starts with the outbound line-haul phase, which goes from the depot to the first client in the district. It is followed by the delivering visits, intercalating cargo discharging tasks with vehicle travel between successive delivery points. Finally, there is the inbound line-haul segment, linking the last served client to the depot.

5.2 Temperature variation modelling

In more elaborated analyses of refrigerated product quality appraisal, the reaction kinetics of the product spoilage process is characterized by determining the kinetics for microbial population growth and chemical degradation. The combined effects of these reactions on perceived organoleptic quality deterioration (appearance, odor, texture, color) is measured and registered in a TTI data sequence (^{Sahin et al., 2007}; Simpson et al., 2011). However, in short-distance distribution cases of refrigerated products, as in the present work, it suffices to register the temperature values in order to verify whether they respect the thermal requisites reflected on PCI values (Section 4).

The theoretical basis to analyze the temperature variation of refrigerated products along time is the Newton's Law of cooling. Suppose that a body, with initial temperature θ ₀, is allowed to cool in a compartment which is maintained at a constant temperature θ ₁, with θ ₁ < θ ₀. Let the temperature of the body be θ at time t . Then, the Newton's law of cooling is

with θ = θ ₀ for t = 0, and where c _s is a positive proportionality constant. Conversely, if θ ₀ < θ ₁, c _s < 0. Suppose a quantity of cargo, with mass m and temperature θ ₀, is set inside an empty refrigerated bin. Departing from the Newton Law of cooling, a simple thermal balance leads to (^{Alvarez & Flick, 1999}; ^{Hoang et al., 2012b})

where C is the specific heat capacity of the product (J kg⁻¹ K⁻¹⁾, H is the heat transfer conductance of the product (W K⁻¹⁾ and β is a constant. The integration of (13) yields an exponential decay function on θ.

In order to fit equations to represent temperature variation along the delivery route, four different phases were considered: (1) outbound line-haul, which corresponds to the vehicle displacement from the depot to the first visited retail outlet; (2) product unloading at retailer premises; (3) vehicle travel from one served retail outlet to the next; (4) inbound line-haul, corresponding to the travel from the last served client back to the depot. A general mathematical expression, to represent the internal vehicle temperature variation along the route, was defined to be fitted to CoolVan simulation results

where

• t _A = starting time of the phase under consideration;

• t _B = time when the phase under consideration ends;

• θ _tA and θ _tB = temperature at times t _A and t _B respectively;

• H = heat transfer conductance of the product, (W K⁻¹⁾;

• C = thermal capacity of the product (J kg⁻¹ k⁻¹⁾;

• m = mass of the product inside the vehicle at the phase under consideration (kg);

• a ₀, a ₁, a ₂ and β are constants to be fitted.

Since in this case the type of product does not vary along all the CoolVan simulation, the values of H and C have been calibrated together with β, and Eq. (14) reduces to

The fitting of constants in Eq. (15) was based on CoolVan simulation values obtained for one vehicle type only, the one described in Table 1. Thus, these results should be understood as preliminary, and should not be extended to other situations. New CoolVan runs are under way, considering other vehicles and other visiting patterns. In the adjusted equations, the mass of the product inside the vehicle in expression (15) is expressed as a fraction of the total vehicle capacity. Furthermore, in the delivery runs analysed in this application, the vehicle leaves the depot with a 100% load, and returns empty to the depot.

a) Outbound line-haul phase (line-haul I): Since the vehicle leaves the depot with a 100% load, m = 1 in Eq. (15). For this phase, considering the vehicle carries the same full cargo load (m = 1) and the external ambient temperature was not statistically significant, Eq. (15) reduces to

Two vehicle displacement intervals were considered: (a) 0 ≤ ΔT ≤ 20 sec, leading to β = −0.7670, and (b) ΔT > 20 sec, leading to β = −0.0044. In addition, a lower limit of 4ºC for the temperature θ _j was imposed for this phase (meaning the refrigeration equipment is turned off when the temperature reaches this level).

b) Cargo unloading phase: a ₀ = 0.3234; a ₁ = 1; a ₂ = 0.4517; β = 0.00081 in relation (15). The product mass varies during the unloading process. Assuming t _B − t _A is the duration of the unloading phase, m _A is the mass at the beginning of the cargo unloading phase, and m _B is the mass at its end, one interpolates linearly over time to determine the value of m. The regression analysis led to R ² = 0.965, and the adjusted constants have shown t Student values significant at a 0.05 level.

c) Vehicle local travel phase: a ₀ = 1.4172, a ₁ = 0.5869; a ₂ = 0.0845 and β = −0.0007 in relation (15). Here m is constant and evaluated at time t _A. The regression analysis led to R ² = 0.913, and the adjusted constants showed t Student values significant at a 0.05 level.

d) Inbound line-haul phase (line haul II): the expression to compute the temperature variation for the inbound line-haul is similar to the outbound line-haul phase. In this phase, the initial temperature is equal to the vehicle air temperature at the moment the last cargo delivery is accomplished. The temperature drops to θ ₁ = 4ºC after an approximately elapsed time of 2 minutes. It is also assumed that no chilled product is sent back to the depot, only working elements such as racks, fittings, empty trays, etc., remain in the truck. An equivalent mass equal to 2.5% of the total load was assumed for such elements, and thus m = 0.025. Applying Eq. (16) for this time interval, one gets β =−0.0112. For T ≥ 2 min to the end of the journey, the temperature remains constant and equal to 4ºC.

5.3 TTI generating module

For a specific routing sequence, the TTI model developed on CoolVan data estimates the travel time between successive delivering stops, the unloading times at different retailer premises, and the line-haul times. Table 2 is an example corresponding to the TSP route. A TTI time series is then constructed assuming one-minute intervals, in which each element is associated with the corresponding phase in the routing process. In the sequel, the vehicle internal temperature is estimated with the equation, indicated in Section 5.2, that corresponds to the phase in which the distribution process is at the moment. An element of the TTI data series is composed by the following elements: (a) element number, (b) phase type, as defined in Section 5.2, (c) starting time t ₀ (min) of the event, (d) termination time t ₁ (min) of the event, (e) temperature θ ₁ (ºC) at time t ₁, estimated with the appropriate equation, (f) last served client number if phase type is local travel, and zero otherwise, (g) next visiting client if phase type is local travel, or the servicing retailer number if it is an unloading phase. Along the computation stages one has for j = 1, 2, . . . and = θ _ini, where θ _ini is the cargo temperature when it is loaded in the vehicle at the depot.

The TTI generationmodulewas first applied to the initial TSP routing sequence shown in Table 2. Since the vehicle returns empty to the depot, the line-haul II phase was not considered in the temperature analysis. A total of 357 TTI pointswere generated. Figure 3 represents the histogram of the temperature inside the vehicle for the TSP routing sequence, showing that its distribution is skewed, which justifies the use of the C _pk index (Section 4).

Figure 3
Histogram of the internal vehicle temperature for the TSP routing alternative.

6 SIMULATED ANNEALING ALGORITHM TO OPTIMIZE THE ROUTING SEQUENCE

6.1 Simulated Annealing

The subject of combinatorial optimization is the development of efficient techniques to find minimum or maximum values of a function with many integer (or binary) variables. Such a function, generally called the objective or cost function, represents a quantitative measure of the performance quality or goodness of a system. Exact methods are computationally unpractical for larger problems with many integer or binary variables and restrictions. In such cases heuristics are usually employed. One basic strategy to solve combinatorial problems is the iterative improvement method. In iterative improvement, one starts with the system in a known configuration, or with an initial solution chosen randomly. A standard rearrangement process is applied to all components of the system until a rearranged configuration, which improves the cost function, is found. The rearranged configuration is now assumed as the new configuration, and the process continues until no further improvements can be found. Iterative improvement comprises a search in the space of feasible solutions for rearrangement steps that will lead downwards toward the global minimum. But this searching process may get stuck in a local minimum. Thus, it is necessary to carry out the process many times, starting from different random configurations and choosing the best solution, which will be taken as the optimum.

In practice, the acceptance of only rearrangements that provoke a reduction in the cost function leads, very often, to poor performance results when applying iterative improvement methods to larger problems. One significant improvement along this line is the Simulated Annealing metaheuristic (^{Kirkpatrick et al., 1983}; ^{Aarts & Korst, 1989}). Simulated annealing (SA) is a generally applicable meta-heuristic for solving combinatorial optimization problems. The algorithmis based on a combination of ideas from two different fields: statistical physics and combinatorial optimization. On the one hand, it can be viewed as an algorithm simulating the physical annealing process of solids to their minimum energy states. On the other hand, it can be considered as a generalization of local search algorithms, which play an important role in combinatorial optimization.

The simulated annealing algorithm starts off with a given initial solution, very often chosen at random, and continuously tries to move from a current solution to one of its neighbours by applying a generation mechanism and an acceptance criterion. The acceptance criterion allows for sporadic deteriorations in the cost function in a limited way. This is controlled by a parameter that plays a similar role as the temperature in the physical annealing process. The possibility of deteriorations makes the SA algorithm more general than pure iterative improvement algorithms, in which only strict improvements are applied. The resulting effect is that the annealing algorithm can systematically avoid local minima in order to arrive at a global minimum (^{Aarts & Korst, 1989}).

^{Metropolis et al. (1953)} developed a simple algorithm that can be used to provide an efficient simulation of a collection of atoms in equilibrium at a given temperature. Using a cost function in place of energy, and defining configurations by a set of parameters {x _i}, ^{Hastings (1970)} generalized the Metropolis procedure to treat iterative improvement optimization problems outside statistical physics, and solving them in an analogous way. Now, temperature is simply a control parameter expressed in the same units as the cost function. In our application, in order to avoid duplication of terms, we call it GT - generalized temperature. The annealing process, as implemented in the Metropolis/Hasting procedure, differs from iterative improvement techniques in that the former avoids getting stuck in local minima, since transitions out of a local optimum are always possible at intermediate temperatures. In general terms, the SA method to solve optimization problems follows the computing sequence depicted in Figure 4.

Figure 4
General SA computing structure.

A number of scientific articles have been published on SA applications to the solution of engineering problems (^{Correia et al., 2000}; ^{Jimenez et al., 2002}; ^{Lin et al., 2009}). Specifically applications on vehicle routing problems, which is the main focus of this article: ^{Aarts, 1989}; ^{Breedam, 1995}; Osman, 1998; ^{Mauri & Lorena, 2006}; ^{Leung et al., 2013}; ^{Hamzadayi et al., 2013}.

6.2 Simulated Annealing applied to the distribution of refrigerated food

Let n be the number of delivering points (clients) to be visited in the route. Each client receives an identification number i = 1, 2, . . . , n. Putting the client numbers according to the sequence they are visited, one has a vector of n components. As mentioned in Section 5, the depot is part of the Hamiltonian cycle, and is placed in position n + 1 on such a vector. This is to permit the direct calculation of the vehicle round trip extension, with data extracted from the distance matrix, including the depot.

Let us initially define a data "type" in pseudo-code:

The data-type contains, first, thementioned array representing the ordering position of the clients, plus the depot. Next it appears the variable C _pk, which is the corresponding PCI value computed for the vector x representing the visiting sequence. Finally, RT is the total route extension, starting from the depot and terminating back to it when all delivering tasks are accomplished. Two variables of type solution are defined: ∅b and ∅a, the first one representing the candidate solution obtained before the SA improvement routine is applied and, the second, the candidate solution generated afterwards.

In the application, two exploratory stages, indicated by k, are considered. In stage k = 1 a candidate solution, represented by the vector x, is defined by a Monte Carlo routine, with the position of each client in the visiting sequence fully generated at random. This is done to allow the SA procedure to explore prospective solutions covering all the set of possibilities, thus avoiding local minima. The sub-optimal solution obtained at the end of stage 1, represented by variable ∅b, will serve as the initial candidate solution for stage 2.

In stage 2, the selection of a new candidate solution is performed in a more restricted manner. In order to get a candidate solution ∅a departing from the former solution ∅b two cients are selected at random in ∅b.x. The positions of these two clients are exchanged mutually in order to generate variable ∅a.x, with the remaining clients keeping the same positions they occupied before in candidate solution ∅b. With this approach the SA improvement sequence will explore closer neighbouring possibilities, thus avoiding too large jumps over the overall set of solutions. Figure 5 shows the pseudo-code of procedure vehicle route (k), that performs the generation of a vehicle route according to the stage k of the SA routine. Other intermediate formulations can also be adopted. For example, one could select a client x _i at random in ∅b.x, and maintain fixed in ∅a.x all the positions of the clients located to the left of x _i in ∅b.x. Next, one shuffles at random all the clients located at the right side of x _i in ∅b.x, and then putting the result in ∅a.x. This scheme brings an intermediate degree of randomness between the other two alternatives. For this application, however, this latter scheme was not necessary because the mathematical structure of the problem led to a rapid SA convergence, as will be explained later.

Figure 5
Pseudo-code of the Monte Carlo procedure to generate vehicle routes.

In applying the SA heuristic to this problem, it is necessary to define the initial value for the generalized temperature (GT ₀ ). Investigating the maximum possible route extension for the problemanalised in this work, one can see that 300 kmis a likely value. With an adequate margin, it was set GT ₀ = 500 km. At each step of the SA meta-heuristic, the generalized temperature is decreased according to a cooling rate rt = 0.99. This selected slow value temperature decreasing rate of 1% allows for a gentle downward trajectory, thus emulating the statistical mechanics annealing process. On the other hand, a minimum generalized temperature level GT _min = 0.01 (practically zero) was adopted.

The initialization part of the SA procedure starts with the definition of an initial vector x _B generated by a Monte Carlo routine, as explained before. With x _B and the distance matrix available for this application, the total route extension RT _B is computed. Similarly, entering the procedure to compute process capability indices (PCI), the corresponding value of is determined (Section 4). The triad {x _B , , RT _B} will constitute variable ∅b.

Let us define a temporary optimum value L∗ for the route extension that is changed whenever a smaller value of RT is encountered. In the initialization process, if < 1.33 (meaning the candidate solution is not acceptable), the value L∗ = GT ₀ is adopted. Otherwise, if ≥ 1.33, L∗ = RT _B is assumed. Next, the computer program enters the SA routine proper. The pseudocode of the Simulated-Annealing procedure is presented in Figure 6.

Figure 6
SA pseudo-code for the refrigerated product distribution problem.

Another candidate solution ∅a is generated using the Monte Carlo routine. If < 1.33, it means that the candidate solution ∅a is not feasible. In this case a penalty is applied duplicating the value of RT _A (row 12 in Figure 6). Otherwise, if ≥ 1.33 and RT _A < L∗, it means that a feasible and better solution has been obtained. In such a case the temporary optimum solution so far is ∅a, and consequently one makes ∅b = ∅a to proceed to a new computational step (rows 15, 16, and 17, Figure 6), and the SA routine reduces the value of the generalized temperature (row 26).

It may occur that the candidate solution ∅a is feasible (with ≤ 1.33) but RT _A ≥ L∗, meaning that the SA routine might have been trapped in a local minimum. In this case, taking a random number 0 ≤ r ≤ 1, one compares the value of r with the analogue Boltzmann probability factor (row 21, in Figure 6). If the inequality in row 21 is observed, the SA accepts ∅a and consequently one makes ∅b = ∅a to proceed to a new computation procedure (rows 22, 23, and 24, Figure 6), and the SA routine reduces the value of the generalized temperature (row 26). Otherwise, the value of the generalized temperature is reduced (row 26) and a new SA iteration is performed. According to this rationale, if the inequality condition in row (21) is satisfied, a worse (i.e., with RT _A ≥ L∗), but feasible solution (i.e., ≥ 1.33) is accepted in order to avoid local minima. Otherwise, the former solution is maintained as provisional optimum so far, and the computation proceeds further to a new iteration.

The stage 1 ends when GT _var ≤ GT _min. In the second stage, the initial candidate solution is not selected at random, but it is just the sub-optimal solution obtained in stage 1.

6.3 SA convergence

As indicated in Section 5.2, four different stage types along the vehicle route are considered: (I) outbound line-haul, which corresponds to the vehicle displacement from the depot to the first visiting retail outlet; (II) product discharge at retailer premises; (III) local vehicle travel from a visited retail outlet to the next one; (IV) inbound line-haul, corresponding to the travel from the last served client back to the depot.

On CoolVan basic results, a mathematical expression to represent the step-by-step internal vehicle temperature variation along the route was calibrated (Section 5.2). In the sequel, a computer program evaluates the vehicle internal temperature at one-minute intervals for each route setting defined in the SA routine. For each route, the simulation produces a complete TTI set, and the resulting C _pk values and total distance D travelled by the vehicle are registered.

The value of constant β in Eq. (15) depends on the type of stage along the route. For stagetypes (I), (III), and (IV) the vehicle is travelling with the cargo door shuttled and the refrigeration system on, and consequently β < 0, since the temperature tends to decrease with time due to the refrigerating process. For stage-type (II) the cargo door stays opened and the refrigeration system is off. The refrigeration system is turned down in such situations to avoid vehicle exhaust pollutants to penetrate the refrigerated compartment, thus spoiling the product. As a consequence β ≥ 0.

One important aspect to be considered in SA applications is the convergence of the iterative process (^{Mengersen & Tweedie, 1996}; ^{Holden, 1998}; ^{Souza de Cursi & Sampaio, 2012}). As mentioned before, in iterative improvement heuristics used to minimize the objective function of a given problem a candidate solution along the process is rejected only when the resulting value of the cost function is not smaller than the previous iterative result. In practice, the dynamics involved in such routines concentrate the search in the neighbourhood of the actual investigated point. The Metropolis-Hastings method, on the other hand, accepts eventual and controlled degradations on the value of the cost function in accordance to a probability evaluation (^{Souza de Cursi & Sampaio, 2012}). From row (21) in Figure 6 one has

^{Souza de Cursi & Sampaio (2012)} presented the following theorem: if {U _i} is a sequence of random variables respecting

such that

the sequence converges to

With decreasing values of GT _var , the probability of accepting a degradation Delta, given by (17), satisfies conditions (18), (19) and (20), leading to the convergence of the Metropolis-Hastings method. In fact, such theorem gives a sufficient but not necessary condition, but the practical result is that theMetropolis-Hastingsmethod has shown to be efficient in solving a large number of practical problems (^{Souza de Cursi & Sampaio, 2012}).

Furthermore, regarding the specific problem treated in this work, the temperature inside the vehicle tends to increase during the product discharge stage, and tends to decrease during the travelling phases. But, since the discharging times are assumed invariable in this analysis, the impact of phase (2) in the optimization process is nil. The same happens with phases (I) and (IV), which also have a fixed behaviour in this analysis. Therefore, the only degrees of freedom when searching for the optimum occurs when the vehicle is travelling from one visit to the next. In other words, the only phase that opens a way to improve the optimization results is phase III (the vehicle travelling phase). Since the mass m of the product inside the vehicle decreases during visiting stops, the value of exp in expression (15) tends to decrease at every SA step, because β ≤ 0 and m has a decreasing value as j increases. In conclusion, the probability of accepting a degradation Delta, given by (17), satisfies conditions (18), (19) and (20), and therefore the SA tends to converge rapidly in our application.

7 RESULTS AND RESEARCH PROSPECTS

7.1 Simulated annealing results

Applying the SA procedures described in Figures 5 and 6, the first SA stage involved 1,306 steps, leading to the following sub-optimal solution:

Routing sequence → Depot - 8 - 10 - 7 - 4 - 3 - 1 - 6 - 12 - 5 - 11 - 9 - 2 - Depot

Total distance travelled → 241.0 km, C _pk = 1.58;

The second SA stage, which departed from the former solution, involved 10,815 steps, leading to the final optimal solution

Routing sequence → Depot - 8 - 7 - 4 - 6 - 9 - 10 - 11 - 3 - 12 - 5 - 2 - 1 - Depot

Total distance travelled → 222.7 km, C _pk = 1.33;

The model was programmed in Free Pascal, version 2.6.4, to be later incorporated in a Delphi format. It runs in a HP Pavilion dm4-2095br notebook. The first and the second stages of the SA algorithm take 4 and 21 sec respectively to run.

The route with minimum vehicle travel distance is the TSP one, with 204.1 km. But this route presented C _pk = 0.35, well bellow the standard level. Taking the inverse TSP route, the PCI is even worse, with C _pk = 0.15. The solution obtained with the SA heuristic has shown an increase of 18.6 km when compared with the TSP configuration, a distance 9.1% greater, but maintaining the refrigerated product within the required temperature range. The SA resulting route is shown in Figure 7, which can be compared with the TSP alternative route. Figure 8, on the other hand, depicts the histogram of the TTI series obtained for the optimal SA solution (as compared with Figure 3, corresponding to the TSP alternative). It can be seen that the temperature distribution is also skewed, justifying the use of the C _pk index (Section 4). The TTI data used to make Figure 8 do not include the values corresponding to the vehicle travel from the last visit back to the depot because the truck is empty along that segment of the route, and therefore the corresponding temperature values do not affect the product quality.

Figure 7
The SA optimized route as compared to the initial TSP route.

Figure 8
SA optimal route: histogram of internal vehicle temperature.

7.2 Sensitivity analysis

For the SA optimal solution indicated in Section 7.1, the values of C _pk were computed for the intermediate stages of the process and shown in Table 3. As mentioned, the stages correspond to the end of the cargo delivery phases. The results of the model show that the thermal properties inside the vehicle are very sensitive to time variations along the route, specifically the cargo discharging times. This happens due to various real problems as, for example, waiting times to start the unloading process, long paths to reach the point where the discharged cargo is to be set inside the client's premises, excessive product checking time, etc. On the other hand, the vehicle travelling times, both along the line-haul segments and along the local links between clients in the route, show comparatively less variation.

Table 3
Optimal SA route: C _pk values for each process stage.

Stage	Client no.	Elapsed time since departure (min)	Cpk
1	8	157	2.16
2	7	185	2.29
3	4	204	2.34
4	6	222	2.39
5	9	239	2.44
6	10	259	2.47
7	11	282	2.46
8	3	306	2.42
9	12	330	2.37
10	5	355	2.23
11	2	376	1.91
12	1	394	1.33

Since there are an expressive number of variables in the problem, it was adopted the coefficient of variation of the vehicle unloading times to perform a preliminary sensitivity analysis of the results, the coefficient of variation (CV ) being the ratio between the standard variation and the expected value. For the outbound and inbound line-haul phases, as well as for the local vehicle displacement times between successive visits, it was assumed CV _LH = 0.1. For the cargo discharging times, on the other hand, a parametric analysis was performed with CV _DT varying from 0.05 to 0.60. Of course, as observed in Table 3, each client shows a particular behaviour in terms of C _pk values. Therefore the sensibility analysis was performed separately for the 12 clients.

Through a Monte Carlo simulation, the model estimates, for each situation, the probability of C _pk < 1.33, as a function of CV _DT , displayed per visiting sequence. In Table 4 and in Figure 9, the probability of not satisfying the thermal requirements are shown. It can be seen in Table 4 and Figure 9 that the evaluation index C _pk values are indeed very sensitive to the product discharging-time variability. The curves exhibit in Figure 8 depend not only on the TTI variation, but also on the quantity of cargo delivered at each point, the unloading time, and other factors. An important conclusion is that the unloading times should be reduced substantially in practice, with extensive field personnel training and the need to perform campaigns directed to retailers, who very often are responsible for unexpected delays in the delivering process. Of course, for a real configuration, with more detailed discharge-time statistical data collected from each client, this analysis can be performed with more accuracy.

Table 4
Probability of occurring C _pk < 1.33 as a function of CV _DT and the visiting sequence.

Visiting sequence	Client number	CVDT - Coefficient of variation, vehicle unloading times
		0.05	0.10	0.15	0.20	0.30	0.40	0.50	0.60
1	8	0	0	0	0	0	0.01	0.02	0.17
2	7	0	0	0	0	0	0.01	0.02	0.17
3	4	0	0	0	0	0	0.02	0.04	0.25
4	6	0	0	0	0	0.01	0.02	0.04	0.25
5	9	0	0	0	0	0.03	0.05	0.09	0.25
6	10	0	0	0	0.01	0.07	0.11	0.17	0.33
7	11	0	0	0	0.03	0.12	0.22	0.27	0.42
8	3	0	0	0.02	0.08	0.21	0.33	0.39	0.58
9	12	0	0.01	0.06	0.13	0.28	0.42	0.50	0.75
10	5	0	0.04	0.09	0.23	0.37	0.49	0.59	0.75
11	2	0.01	0.11	0.16	0.37	0.49	0.61	0.67	0.75
12	1	0.31	0.40	0.49	0.54	0.62	0.72	0.75	0.83

Figure 9
Probability of not satisfying the thermal requirement C _pk ≥ 1.33, as a function of CV _DT .

7.3 Conclusion and research prospects

The main contribution of this work is the joint investigation of the thermal properties of a refrigerated product along a distribution route, where TTI data sets are gathered and the product quality is evaluated with a process capability index. Alternative routing sequences are then optimized with a Simulated Annealing algorithm in which the objective function to be minimized is the vehicle travelled distance, but respecting a minimum capability index target value.

The authors intend to research further on this subject considering a larger refrigerated cargo distribution problem and using an irregular fleet of vehicles, including smaller trucks, with each one performing more than one route in a working day, the vehicles getting back to the depot when a visiting tour is finished and performing other visiting delivery sequences, as in ^{Oswald & Stirn (2008)}. Such scheme will possibly improve the overall thermal performance of the delivery system, although with somewhat higher operating costs.

Second, additional TTI data will be gathered with the CoolVan simulator, in order to improve the regression fittings representing the TTI variation. Finally, a dynamic version of this coldchain distribution planning problem will be developed considering other possible time deviations that might deserve routing revisions. For example, unexpected delays in the cargo unloading activities, traffic delays, and other fault occurrences that will be dynamically incorporated into the model in such a way as to devise corrective real-time measures.

ACKNOWLEDGMENTS

This research has been supported by Capes Foundation (Brazil) and DFG, the German Research Foundation, Bragecrim Project #2009-2. We sincerely thank one of the reviewers for the thorough and helpful revision of the paper.

References

• ¹AARTS EHL & KORST JHM. 1989. Boltzmann machines for travelling salesman problems. European Journal of Operational Research, 39: 79-95.
• ²ALVAREZ G & FLICK D. 1999. Analysis of heterogeneous cooling of agricultural products inside beans Part II: thermal study. Journal of Food Engineering, 39: 239-245.
• ³ANDRÉE S, JIRA W, SCHWIND K-H, WAGNER H & SCHWÄGELE F. 2010. Chemical safety of meat and meat products.Meat Science, 86: 38-48.
• ⁴BORCH E & ARINDER P. 2002. Bacteriological safety issues in red meat and ready-to-eat meat products, as well as control measures.Meat Science, 62: 381-390.
• ⁵BREEDAM AV. 1995. Improvement heuristics for the vehicle routing problem based on simulated annealing. European Journal of Operational Research, 86: 480-490.
• ⁶CAMPAÑONE LA, GINER SA & MASCHERONI RH. 2002. Generalized model for the simulation of food refrigeration. Development and validation of the predictive numerical method. International Journal of Refrigeration, 25: 975-984.
• ⁷CHANG YS & BAI DS. 2001. Control charts for positively-skewed populations with weighted standard deviations. Quality and Reliability Engineering International, 17: 397-406.
• ⁸CHANG YS, CHOI IS & BAI DS. 2002. Process capability indices for skewed populations. Quality and Reliability Engineering International, 18: 383-393.
• ⁹CHEN JP & DING CG. 2001. A new process capability index for non-normal distributions. International Journal of Production Research, 38(6): 1311-1324.
• ¹⁰COOLVAN MANUAL - VERSION 3.0. 2000. FRPERC - Food Refrigeration & Process Engineering Research Centre, University of Bristol, UK.
• ¹¹CORREIA MH, OLIVEIRA JF & FERREIRA JS. 2000. Cylinder packing by simulated annealing. Pesquisa Operacional, 20(2): 269-286.
• ¹²CUESTA JF, LAMÚA M & MORENO J. 1990. Graphical calculation of half-cooling times. Int. J. Refrigeration, 13: 317-324.
• ¹³CZARSKI A. 2008. Estimation of process capability indices in case of distribution unlike the normal one. Archives of Materials Science and Engineering, 34(1), 39-42.
• ¹⁴DAELMAN J, JACXSENS L, DEVLIEGHERE F & UYTTENDAELE M. 2013. Microbial safety and quality of various types of cooked chilled foods. Food Control, 30: 510-517.
• ¹⁵ELLOUZE M & AUGUSTIN J-C. 2010. Applicability of biological time temperature integrators as quality and safety indicators for meat products. International Journal of Food Microbiology, 138: 119-129.
• ¹⁶ESTRADA-FLORESS & TANNER D. 2005. Interpretation of the Australian Standard for the thermal testing of refrigerated vehicles. EcoLibrium, November 2005, pp. 24-31.
• ¹⁷ESTRADA-FLORESS & EDDY A. 2006. Thermal performance indicators for refrigerated road vehicles. Int. J. Refrigeration, 29: 889-898.
• ¹⁸FARAZ A, HEUCHENNE ES & FOSTER E. 2013.Monitoring delivery chains using multivariate control charts. European Journal of Operational Research, 228(1): 282-289.
• ¹⁹FERREIRA VO & PUREZA V 2012. Some experiments with a savings heuristic and a tabu search approach for the vehicle routing problemwith multiple deliveries. Pesquisa Operacional, 32(2): 443-463.
• ²⁰FLICK D, HOANG HM, ALVAREZG & LAGUERREO. 2012. Combined deterministic and stochastic approaches for modelling the evolution of food products along the cold chain. Part I: Methodology. Int. J. Refrigeration, 35: 907-914.
• ²¹GARCIA LR. 2008. Development of Monitoring applications for refrigerated perishable goods transportation. PhD Thesis, Universidad Politécnica de Madrid, Spain.
• ²²GE YT & TASSOU SA. 2001. Simulation of the performance of single jet air curtains for vertical refrigerated display cabinets. Applied Thermal Engineering, 21, pp. 201-219.
• ²³GIANNAKOUROU MC & TAOUKIS PS. 2003. Application of a TTI-based distribution management system for quality optimization of frozen vegetables at the consumer end. Journal of Food Science, 68(1): 201-209.
• ²⁴GIANNAKOUROU MC, KOUTSOUMANIS K, NYCHAS GJ & TAOUKIS PS. 2005. Field Evaluation of the application of time temperature integrators for monitoring fish quality in the chill chain. Int. J. Food Microbiology, 102: 323-336.
• ²⁵GIGIEL A. 1997. Predicting food temperature in refrigerated transport. Proceedings: The Institute of Refrigeration, UK: 7-1 to 7-12.
• ²⁶GIGIEL A, JAMES SJ & EVANS JA. 1998. Controlling temperature during distribution and retail. Proceedings 3rd Karlsruhe Nutrition Symposium: European Research towards Safer and Better Food, Karlsruhe, Germany, 284-292.
• ²⁷GONÇALEZ PU & WERNER L. 2009. Comparison of process capability indices for non-normal distributions (in Portuguese). Gestão & Produção, 161: 121-132.
• ²⁸HAMZADAYI A, TOPALOGLUS & KOSE SY. 2013.Nested simulated annealing approach to periodic routing problem of a retail distribution system. Computers & Operations Research, 40(12): 2893-2905.
• ²⁹HARVEY JM. 1978. Reduction of losses in fresh market fruits and vegetables. Ann. Rev. Phytopathology, 16: 321-341.
• ³⁰HASTINGS WK. 1970. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57: 97-109.
• ³¹HAVELAAR AH, BRUL S, DE JONG A, DE JONGE R, ZWIETERING MH & TER KUILE B. 2010. Future challenges to microbial food safety. International Journal of Food Microbiology, 139: S79-S94.
• ³²HOANG MH, FLICK D, DERENS E, ALVAREZ G& LAGUERRE O. 2012a. Combined deterministic and stochastic approaches for modelling the evolution of food products along the cold chain. Part II: A case study. Int. J. Refrigeration, 35: 915-926.
• ³³HOANG MH, LAGUERRE O, MOURETH J & FLICK D. 2012b. Heat transfer modelling in a ventilated cavity loaded with food product: Application to a refrigerated vehicle. J. of Food Engineering, 113: 389-398.
• ³⁴HOLDEN L. 1998. Geometric convergence of the Metropolis-Hastings simulation algorithm, Norwegian Computing Center, P.O. Box 114, Blindern, N-0314, Oslo, Norway.
• ³⁵ISO 1496-2 E, SERIES 1. 1996. Freight Containers - Specification and Testing - Part 2: Thermal Containers.
• ³⁶JACXSENS L, DEVLIEGHERE F& DEBEVERE J. 2002. Predictive modelling for packaging design: equilibrium modified atmosphere packages of fresh-cut vegetables subjected to a simulated distribution chain. Int. J. Food Microbiology, 73: 331-341.
• ³⁷JAMES SJ & SCHOLFIELD I. 1998. Modelling of food refrigeration systems. Proceedings 3rd Karlsruhe Nutrition Symposium: European Research towards Safer and Better Food,Karlsruhe, Germany, 293-301.
• ³⁸JAMES SJ, JAMES C & EVANS JA. 2006. Modelling of food transportation systems - a review. Int. J. Refrigeration, 29: 947-957.
• ³⁹JAMES SJ & JAMES C. 2010. The cold chain and climate change. Food Research International, 43: 1944-1956.
• ⁴⁰JIMENEZ A, RIOS-INSUAS & MATEOS A. 2002. Interactive simulated annealing for solving imprecise discrete multiattribute problems under risk. Pesquisa Operacional, 22(2): 265-280.
• ⁴¹KIRKPATRICK S, GELLAT CD & VECCHI MP. 1983. Optimization by simulated annealing. Science, 220: 671-680.
• ⁴²LEUNG SCH, ZHANG Z, ZHANG D, HUA X & LIM MK. 2013. A meta-heuristic algorithm for heterogeneous fleet routing problems with two-dimensional loading constraints. European Journal of Operational Research, 225: 199-210.
• ⁴³LIN S-W, YU VF & CHOU S-Y. 2009. Solving the truck and trailer routing problem based on a simulated annealing heuristic. Computers & Operations Research, 36: 1683-1602.
• ⁴⁴MAI NT, GUDJÓNSDÓTTIR M, LAUZON HL, SVEINSDÓTTIR K, MARTINSDÓTTIR E, AUDORFF H, REICHSTEIN W, HAARER D, BOGASON SG & ARASON S. 2011. Continuous quality and shelf life monitoring of retail-packed fresh cold loins in comparison with conventional methods. Food Control, 22: 1000-1007.
• ⁴⁵MAURI GR & LORENA LAN. 2006. Simulated annealing applied to a general model to solve the vehicle routing problem (in Portuguese), XXXVIII SBPO - Brazilian Symposium of Operations Research, Goiânia, GO, September, 12 to 15, 2006, Brazil.
• ⁴⁶MCMEEKIN T, BOWMAN J, MCQUESTIN O, MELLEFONT L, ROSS T & TAMPLIN M. 2008. The future of predictive microbiology: Strategic research, innovative applications and great expectations. International Journal of Food Microbiology, 128: 2-9.
• ⁴⁷MENDOZA TF, WELT BA, OTWELL S, TEIXEIRA AA, KRISTONSSONH & BALABAN MO. 2004. Kinetic parameter estimation of time-temperature integrators intended for use with packaged fresh seafood. Journal of Food Science, 69(3): FMS90-FMS96.
• ⁴⁸MENGERSEN KL & TWEEDIE RL. 1996. Rates of convergence of the Hastings and Metropolis algorithms. The Annals of Statistics, 24(1): 101-121.
• ⁴⁹METROPOLIS N, ROSENBLUTH A, ROSENBLUTH M, TELLER A & TELLER E. 1953. Equations of state calculations by fast computing machines. Journal of Chemical Physics, 21: 1087-1091.
• ⁵⁰MOURETH J & DERENS E. 2000. Numerical modelling of the temperature increase in frozen food packaged in pallets. Int. J. Refrigeration, 23: 540-552.
• ⁵¹MOUREH J, MENIA N & FLICK D. 2002. Numerical and experimental study of airflow in a typical refrigerated truck configuration loaded with pallets. Computers and Electronics in Agriculture, 34: 25-42.
• ⁵²NBR 14701. 2001. Transport of refrigerated foodstuffs - Procedure and criterion of temperature, ABNT - Brazilian Association of Technical Norms, Brazil.
• ⁵³NGA MT. 2010. Enhancing Quality Management of Fresh Fish Supply Chains Through Improved Logistics and Ensured Traceability, PhD Thesis, University of Iceland, Reykjavik.
• ⁵⁴OSMAN IH. 1993. Metastrategy simulated annealing and tabu search algorithms for the vehicle routing problem. Annals of Operations Research, 41: 421-451.
• ⁵⁵OSWALD A & STIRN LZ. 2008. A vehicle routing algorithm for the distribution of flesh vegetables and similar perishable food. Journal of Food Engineering, 85: 285-295.
• ⁵⁶PEARN WL & CHEN KS. 1999.Making decisions in Assessing Process Capability Index Cpk. Quality and Reliability Engineering International, 15: 321-326.
• ⁵⁷PEREIRA VF, DÓRIA EC, CARVALHO JR BC, NEVES JR LC, & SILVEIRA JR V. 2010. Evaluation of temperatures in a refrigerated container for chilled and frozen food transport (in Portuguese). Ciência e Tecnologia de Alimentos, Campinas, 30(1), 158-165.
• ⁵⁸SAHIN E, BABAÏ MZ, DALLERY Y & VAILLANT R. 2007. Ensuring supply chain safety through time temperature integrators. Int. J. Logistics Management, 18(1): 102-124.
• ⁵⁹SIMPSON R, ALMONACID S, NUÑEZ H, PINTO M, ABAKAROVA & TEIXEIRA A. 2012. Time-temperature indicator tomonitor cold chain distribution of fresh salmon (salmo salar). J. Food Process Engineering, 35(5): 742-750.
• ⁶⁰SMOLANDER M, ALAKOMI H-L, RITVANEN T, VAINIONPÄÄ J & AHVENAIMER R. 2004. Monitoring of the quality of modified atmosphere packaged broiler chicken cuts stored in different temperature conditions. A. Time-temperature indicators as quality-indicating tools. Food Control, 15: 217-229.
• ⁶¹SOUZA DE CURSI JE & SAMPAIO R. 2012. Stochastic modelling and uncertainties quantification (in Portuguese). Brazilian Association of Applied and Computational Mathematics, São Carlos, SP, Brazil.
• ⁶²SYSLO MM, DEO N & KOWALIK JS. 2006. Discrete Optimization Algorithms with Pascal Programs. Dover Publications, Mineola, NY.
• ⁶³TARANTILIS CD & KIRANOUDIS CT. 2002. Distribution of fresh meat. Journal of Food Engineering, 51: 85-91.
• ⁶⁴TASSOU SA, DE-LILLE G & GE YT. 2009. Food transport refrigeration - Approaches to reduce energy consumption and environmental impacts of road transport. Applied Thermal Engineering, 29: 1467-1477.
• ⁶⁵TASSOU AS, DE-LILLE G& LEWIS J. 2012. Food transport refrigeration. Centre for Energy and Built Environment Research, School of Engineering and Design, Brunel University, UK.
• ⁶⁶TSIRONI T, STAMATIOU A, GIANNOGLOU M, VELLION E & TAOUKIS PS. 2001. Predictive modelling and selection of time temperature integrators for monitoring the shelf life of modified atmosphere packed gilthead seabreamfillets. LWT - Food Science and Technology, 44: 1156-1163.
• ⁶⁷TSO CP, YU SC, POH HJ & JOLLY PG. 2002. Experimental study of the heat and mass transfer characteristics in a refrigerated truck. Int. J. Refrigeration, 25: 340-350.
• ⁶⁸WEGENER HC. 2010. Danish initiatives to improve the safety of meat products. Meat Science, 84: 276-283.
• ⁶⁹YU M & NAGURNEY A. 2013. Competitive food supply chain networks with application to fresh produce. European Journal of Operational Research, 224: 273-282.
• ⁷⁰ZHANG G & SUN G. 1994. A new method to determine the heat transfer coefficient of refrigerated vehicles. Int. J. Refrigeration, 17(8): 516-523.
• ⁷¹ZHOU GH, XU XL & LIU Y. 2010. Preservation Technologies of fresh meat - A review. Meat Science, 86: 119-128.

Publication Dates

• Publication in this collection
  May-Aug 2015

History

• Received
  27 Nov 2013
• Accepted
  12 Aug 2014