Process capability index Cpk for monitoring the thermal performance in the distribution of refrigerated products

The temperature of refrigerated products along the cold chain must be kept within pre-defined limits to ensure adequate safety levels and high product quality. Because temperature largely influences microbial activities, the continuous monitoring of the time-temperature history over the distribution process usually allows for the adequate control of the product quality along both shortand medium-distance distribution routes. Time-Temperature Indicators (TTI) are composed of temperature measurements taken at various time intervals and are used to feed analytic models that monitor the impacts of temperature on product quality. Process Capability Indices (PCI), however, are calculated using TTI series to evaluate whether the thermal characteristics of the process are within the specified range. In this application, a refrigerated food delivery route is investigated using a simulated annealing algorithm that considers alternative delivery schemes. The objective of this investigation is to minimize the distance traveled while maintaining the vehicle temperature within the prescribed capability level.


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 similar products.In order to ensure product quality and health safety, 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 the temperature at which the product is held during handling, storage, transport, and distribution, the use of controlled or modified atmospheres during storage or transit, chemical treatments for the control of decay or physiological disorders, heat treatments for decay control, packaging and handling systems, etc. (James & James, 2010).But, since temperature largely determines the rate of microbial activity, which is the main cause of spoilage of most fresh food products (Borch & Arinder, 2002), continuous monitoring of the full time temperature history usually allows for an adequate control of the process along the short and medium distance distribution routes (Simpson et al., 2012;Flick et al., 2012).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.Metabolic activities generally increase as storing temperatures are elevated.On the other hand, short interruptions in the control of the cold chain may result quick deterioration of product quality (Jedermann et al., 2009).Hence, the required temperature range should be maintained from production to consumption.
In addition, along 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 (James et al., 2006).
A number of papers on refrigerated food transport address TTI -time-temperature evaluation along the cold chain (Gigiel et al., 1998;Giannakourou & Taoukis, 2003;Smolander et al., 2004;Giannakourou et al., 2005;Estrada-Flores & Eddy, 2006;Sahin et al., 2007;Pereira et al., 2010).For short and medium distance product delivery, 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, in fact, adopted in a number of studies reported in the literature and based on TTI.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.
Process Capability Indices (PCI), on the other hand, can yield easily computed coefficients measured with dimensionless functions on TTI parameters and specifications.(Chang et al., 2002;Barriga et al., 2003;Bulba & Ho, 2004;Mingoti & Glória, 2008;Chang, 2009;Gonçalez & Werner, 2009, Mingoti et al., 2011).This kind of PCI application helps to reveal undesired thermal conditions that may impair the compliance of product quality requirements along the supply chain.
This paper reports a Time-Temperature Indicator (TTI) analysis of a distribution of refrigerated food products along a route containing a number of retail customers with different demand levels.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 is minimized, usually leads to temperature restriction violations.Thus, instead of using a classical method to get the optimized TSP vehicle routing sequence, such as the largely employed 2-opt and 3-opt improvement heuristics (Syslo et al., 2006), other solving algorithms, apart from the TSP, can be employed as, for example, the Simulated Annealing (SA) meta-heuristic (Breedam, 1995;Mauri & Lorena, 2009;Leung et al., 2013).In this paper, a SA algorithm was specifically developed to optimize the routing problem.

Related works and research contribution
A number of papers on refrigerated food transport address time-temperature evaluation along the cold chain process (Gigiel et al., 1998;Giannakourou et al., 2005;Estrada-Flores & Eddy, 2006;Pereira et al., 2010;Simpson et al., 2012).TTI data for such studies are usually gathered in field surveys (Pereira et al., 2010), in performance laboratory tests (Moureh & Derens, 2000;Tso et al., 2002;Estrada-Flores & Eddy, 2006;Jedermann et al., 2009), or with the help of computer based simulations (Gigiel et al., 1998;Food Refrigeration & Process Engineering Research Centre, 2000).Considering the objectives of this study, which involves a practical logistics application, the thermal laboratory testing alternative 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 process, which only requires a simpler data approximation.Field tests to gather TTI data, on the other hand, could only be considered under rigid technological and operational control, a condition hard to achieve presently in developing countries like Brazil (Pereira et al., 2010).Thus, the third alternative was chosen for this work, and the CoolVan simulation software was used to gather TTI data (Section 3).
Vehicle routing problems in the distribution of perishable food were previously analysed by Tarantilis & Kiranoudis (2002), Hsu et al. (2007), Azi et al. (2007), Oswald & Stirn (2008), Chen et al. (2009), andAmorim et al. (2012).Tarantilis & Kiranoudis (2002) developed an algorithm to optimize the distribution of fresh meat formulated as a multidepot vehicle routing problem carrying the product to butchers' shops located in an area of the city of Athens, Greece.Product temperature variation along the delivering process was not explicitly analysed in the paper, implying that thermal behaviour was not critical under the conditions involved in their work.Thus, Tarantilis & Kiranoudis (2002) paper can be regarded as a typical application of a classical multidepot vehicle routing problem (VRP).route, searching for the minimum travel distance configuration, but at same time respecting a lower bound value for the chosen quality index.

TTI data obtained via CoolVan simulation
One of the most systematic attempts to predict the temperature of refrigerated food during multidrop deliveries has been the CoolVan research programme.The CoolVan software was developed by the Food Refrigeration and Process Engineering Research Centre, at the University of Bristol, UK (Gigiel et al., 1998;Food Refrigeration & Process Engineering Research Centre, 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 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.
The usual food distribution scheme starts with the vehicle being loaded at the distributor's depot and travelling 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.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 for 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.Presently, the CoolVan software is not available commercially, but its developers kindly Hsu et al. (2007), Oswald & Stirn (2008), Chen et al. (2009), andAmorim et al. (2012) regarded the vehicle routing problem of perishable products with a different quality assurance criterion.Instead of TTI control, as in this paper, these authors considered a time-sensitive spoilage rate of the product.Particularly Hsu et al. (2007) defined an interesting spoilage evaluation procedure in which, as time goes by along the vehicle visiting sequence, the estimated fraction of spoiled product increases.It is assumed that the spoilage of food has not yet begun at the distribution centre, when the product has been loaded into the vehicle.The spoilage rate increases at a certain level when the truck door is opened for discharge.Spoilage is also observed at a different rate when the vehicle travels from one customer to the next.Oswald & Stirn (2008) adopted a similar, but simpler quality criterion, in which product quality (fresh vegetables) is based on market acceptance: one has 100% quality when 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.
In the literature, no work was encountered in which TTI data were used in association with a mathematical vehicle routing model involving perishable products.A research contribution of this work has been the integration, into a same modelling framework, the TTI data analysis, associated with a process capability index evaluation, and using a SA meta-heuristic to optimize the vehicle routing sequence, but keeping the thermal performance indicators within the required levels.In addition, a computer framework analysis was specifically developed to solve the problem, involving a combination of methods divided in four steps: • First, the CoolVan software was used to simulate the thermal characteristics of the basic routing schemes, one of them representing the TSP formulation.
• Next, taking the CoolVan simulation results, the operating stages that compose the TTI evolution along the distribution journey (line-haul linking the depot to the district and vice-versa, cargo discharging at the retailers' premises, and local vehicle travel) were analysed individually in order to define mathematical functions that relate temperature to the explaining variables.
• Third, a computer program, specially developed to estimate TTI values step by step along each routing sequence, calculates the corresponding process capability index.This index allows the comparison of the thermal performance of each alternative delivery sequence.
• Finally, a simulated annealing algorithm (Section 5) was developed to analyse each candidate 4.1.Capability indices C p , C pk , C pm and C pmk 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 distribution journey of refrigerated food products, with mean µ and standard deviation σ.These are unknown elements and are estimated by µ ^ and σ ^ respectively.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.Table 1 indicates quality status as a function of C p values, according to Kaya & Kahraman (2010).
The six-sigma coverage represents the spread of 99.73% of the data in normally distributed processes.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 (2) The quality conditions indicated in Table 1 for C p are also valid for C pk values.On the other hand, when the distribution of θ is normal and the mean offered to ran some configurations to serve as a data basis for this application.

Process capability indices in thermal performance analysis
A process capability index is a numerical element that compares the characteristics of a production or servicing process to engineering specifications.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.A capability index of this sort is convenient because it reduces complex information about the quality of the process to a single number (Pearn & Lin, 2004;Anis, 2008;Wu et al., 2009;Gonçalez & Werner, 2009).A former application of the PCI methodology in the thermal evaluation of the transport of refrigerated products is described in Estrada-Flores & Eddy (2006).
The starting point of the process capability analysis is the definition of a few measurable properties which can give a significant insight into the quality of the output of the process under consideration.For that, one has to associate a specification interval, that is, an upper specification limit (USL) and a lower specification limit (LSL), to each of the measurements of interest.Of course, some measurement values might be out of the specification limits.Thus, the capability analysis problem consists in the estimation of a statistical model for the process in order to be able to predict the number of points falling out of the specification limits, and therefore giving a measure of the process capability (Bittanti et al., 1998).
Process capability indices (PCI) are frequently used as an integral part of the statistical control of system quality and productivity.The relationship between the actual process performance and the specification limits or tolerance may be quantified using appropriate PCI.They are statistically designed to provide common and easily computed coefficients measured with a dimensionless function of its parameters and specifications.In our application, the analysis of the temperature variability inside a refrigerated vehicle is performed by means of TTI data series collected along typical distribution journeys.Since temperature largely determines the rate of microbial activity, continuous monitoring of the full temperature history usually allows for the adequate control of the process along short and medium distance situations.families of frequency curves to select an adequate probability distribution to that objective.The method of moments is a known approach to the problem of estimating the parameters from a family of frequency curves (Bittanti et al., 1998;Elderton & Johnson, 1969).It is based on the following: given a family of probability distributions which depends on a number of parameters, one can manage to derive close form expressions for the first moments, by applying the operators to the probability density function f (θ, ϕ) of θ, where ϕ ∈ R n is the vector of unknown parameters.Then, by equating such parametric expressions to the numerical estimates of the moments obtained from the data, one has a set of n equations which can be solved for the n unknown parameters.For common families of distributions, the number of parameters is not very large.Thus, it is usually possible to solve the equations analytically and obtain transformations with the estimated moments as inputs, and the values of the parameters as outputs (Bittanti et al., 1998).
The Pearson system of curves is defined as the set of functions f (θ) satisfying the following differential equation (Bittanti et al., 1998;Elderton & Johnson, 1969): where a, b 0, b 1 and b 2 are four real parameters.Among the solutions of ( 7), curves of different kinds can be found, like, for example, bell-shaped, J-shaped or U-shaped distributions.The solution of (7) varies according to the roots of the second order polynomial (Elderton & Johnson, 1969): As proposed by Pearson, the various forms or curve types can be classified according to the value of the so-called critical element K: which ranges from -∞ to +∞ and is related to the location of the roots of polynomial (8) in the complex plane.
Another form to express K is (Bittanti et al., 1998) ( ) where µ ^ is centred on the target value T, one has C p = C pk .
In addition, the index C pm explicitly considers the distance between the mean and the target value T (Gonçalez & Werner, 2009): To obtain a capability index more sensitive than C pk and C pm , with regard to departures of the process mean µ ^ from the target value T, another third generation of capability index C pmk has been introduced (Gonçalez & Werner, 2009;Chang, 2009): 4.2.Capability analysis of non-normal data

Clements' method
When it is detected that the process control variable is not normally distributed, one way to introduce the effects of non-normality into the analysis is Clements' method, in which the expression (1) is generalized as follows (Bittanti et al., 1998;Gonçalez & Werner, 2009;Hosseinifard et al., 2009): where m is the median of the distribution and θ 0.9987 and θ 0.0013 are the percentiles corresponding to probabilities 0.9987 and 0.0013, respectively.It is clear that in the normal case this formulation reduces to (2), since then µ ^ = m, θ 0,0013 = µ ^ -3σ ^, and θ 0,9987 = µ ^ + 3σ ^.This nonparametric method seems very attractive because of its simplicity, but unfortunately it suffers from two major drawbacks.First, when working with small samples or noisy data, the nonparametric estimation of tail percentiles becomes a difficult task and the obtained estimates can be very inaccurate.Therefore, this approach is likely to be very sensitive to sampling in actual applications.Second, the estimated percentiles do not provide a very clear picture of the distribution underlying the data.This signifies that it cannot serve for process monitoring, because of its limited ability to detect changes in the process distribution (Bittanti et al., 1998).

The Pearson system of curves
Another approach to handle non-normal situations is to search for a type of probability distribution that best fits to the data.One way to do this is to base the analysis on parametric In our application, the variable θ is not represented by a normal distribution (see Section 6).As indicated in Section 6, the TTI data obtained with the SA meta-heuristic are represented by 394 temperature and time values covering a typical tour, to the exception of the line-haul segment from the last visit to the depot.Computing the central moments of the TTI series, one has m 2 = 2.227, m 3 = 0.352 and m 4 = 22.922, with mean µ = 5.05 °C and skewness = 0.106.Expression (10), associated with (11), yields K = 0.0062 > 0, meaning that the probability distribution of θ is a Pearson's Type IV curve.The fitting of a distribution such as ( 14) is not simple and usually requires large samples in order to adequately represent the tail effect in the mathematical representation.

The Burr transformation approach
Another method of handling non-normal situations in PCI analysis is to reduce the non-normal data into a normal configuration through a transformation technique, and then use a conventional normal method to estimate PCI values (Lin & Chen, 2006).Burr (1973) proposed a system of probability distributions, denominated Burr XII distribution, whose cumulative function is The corresponding probability density function of the Burr XII distribution is (Lin & Chen, 2006) (17) Burr (1973) tabulated the expected values, standard deviations, skewness and kurtosis coefficients for the Burr XII distribution, for several combinations of c and k.The resulting tables enable users to make a standardized transformation between a sample variate and a Burr random variate.For a set of data, after the sample skewness and kurtosis coefficients have been estimated, the mean and standard deviation of the corresponding Burr distribution may be obtained using standard tables.
As mentioned, the fitting of probability functions, such as Pearson's curves and Burr's distribution, is not simple and usually requires large samples in order to adequately represent the tail effect in the mathematical representation.As a result, several approximate approaches to the PCI problem with skewed populations have been proposed and tested in the literature (Chang & Bai, 2001;Chang et al., 2002;Gonçalez & Werner, 2009).Pearson defined 13 such types of curves but, for practical applications, three main types are of interest (Bittanti et al., 1998;Elderton & Johnson, 1969): Pearson's Type I: Expression (8) has two real roots of opposite sign, yielding K < 0. The solution of ( 8) is ( ) where d 1 , d 2 , g 1 , g 2 are coefficients to be fitted to the data and y 0 is the normalization factor.The corresponding pdf of ( 12) only exists when θ is limited to the range over which θ is real and positive.
Pearson's Type IV: It corresponds to values of K between 0 and 1, i.e. to situations where (8) has two complex conjugate roots, leading to the solution ( ) where d, ν, r, g are the four parameters to be fitted to the data and y 0 is the normalization factor.The range of θ in this case is unlimited, i.e. -∞ < θ <+∞.
Pearson's Type VI: It corresponds to values of K greater than one, i.e., to situations where (8) has two real roots of the same sign, leading to the solution ( ) with the pdf of f (θ) limited to the range d 1 < θ <+∞, or -∞ < θ < d 2 , depending on the sign of the roots of polynomial (8).Thus, while Type I curves are adequate for the representation of variates with both an upper and a lower limit, Type VI curves can be used when only one such limit exists as, for example, the case of a pdf of a nonnegative variable to be fitted to the data (Bittanti et al., 1998).All other Pearson's curve types, sometimes indicated as transition types, correspond to single values of K.For instance, when K is nil, it indicates that the distribution is symmetrical and, in the Pearson system of frequency curves, it can be the Normal, the Uniform, the Pearson type II and the Pearson type VII.
The PCI for temperature θ, for the case of two-side specifications (upper and lower values) is In Equations 22 and 23, 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 ( 24) is smaller than the standard one given by (1).
The numerical estimation of C p and C pk is performed as follows.Let (θ 1 , θ 2 , ..., θ 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 with value smaller 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 P ^θ for µ, σ and P θ respectively in Equation 24.

A simulated annealing approach to solve the vehicle routing problem
Simulated annealing (SA) is a meta-heuristic for solving combinatorial optimization problems.The algorithm is based on a combination of ideas from two different fields: statistical physics (physical annealing process of solids converging to their minimum energy states) and combinatorial optimization (Kirkpatrick et al., 1983;Magalhães & Moura Neto, 2011).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 optimization 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.
Many works have been published on SA applications to the solution of vehicle routing 4.3.An approximate C pk evaluation method 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 better accuracy the number of normal non-conforming items in the evaluated sample, being superior to the other analysed 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 ( ) ( ) 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) ( ) with P θ = Pr{θ < µ}.
The value of C p is defined as (Chang & Bai, 2001) 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 TTI data used to make Figure 2 do not include the values corresponding to the vehicle travel from the last visit back to the depot because, in this application, the truck travels empty along the backwards line-haul segment of the route since all cargo has been delivered up to that point, and therefore the corresponding temperature values do not affect the product quality.The mean temperature is µ ^ = 5.05 °C, and its standard deviation σ ^ = 1.49.problems (Breedam, 1995;Osman, 1993;Mauri & Lorena, 2009;Leung et al., 2013).In this paper, SA has been applied to solve the optimization of the proposed vehicle routing problem.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.The depot is part of the Hamiltonian routing cycle, and is placed in position n + 1 on such a vector.In the application, two exploratory stages, represented 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 selected at random.This is done to allow the SA procedure to explore prospective solutions covering all the set of possibilities, thus avoiding local minima.In stage 2, the selection of a new candidate solution is performed in a more restricted manner.In order to get a candidate solution departing from a former solution, two clients are selected at random in the corresponding vector.The positions of these two clients are exchanged mutually in order to generate a new prospective solution, with the remaining clients keeping the same positions they occupied before.With this approach the SA improvement sequence explores closer neighbouring possibilities, thus avoiding too large jumps over the overall set of solutions.

Results and sensitivity analysis
The necessary basic inputs for the model application are presented in Table 2.The urban distribution district in this work 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.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 route starts with a line-haul phase, which goes from the depot to the first client in the district, and taking 84 minutes approximately.It is followed by a sequence of visits, intercalating cargo discharging tasks with vehicle displacements between successive delivery points.Finally, there is the line-haul reverse segment, linking the last served client to the depot.
The optimal delivery sequence shown in Figure 1 was obtained with the SA algorithm, in which the objective is to minimize the vehicle travelling  is even worse, with C pk = 0.15.Looking at the zigzag sequence pattern of the SA solution (Figure 1), it might look non intuitive, but the increased distances linking some servicing points can help the vehicle refrigeration equipment to reduce the internal temperature to acceptable levels.In fact, the temperature inside the truck along the delivery journey tends to increase when it stops to discharge the product, and the temperature inversely tends to decrease while the vehicle is travelling from a client to the next one in the routing sequence.Generally, a heavier discharge weight takes more time to be unloaded.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 in the TTI history, indicating that both have to be investigated jointly, as it has been done in this analysis.
For the SA optimal solution indicated above, the values of C pk were also computed considering the intermediate stages of the process.The process stages correspond to the instants when the vehicle has just finished a delivery service to a client, and is preparing to depart for the next visit.In order to investigate the sensitivity of these results, a simulation was performed considering the product discharging times and the vehicle travelling times as random variables.For this, such variables were assumed to be log-normally distributed.It was assumed a coefficient of variation CV = 0.2 for the cargo discharging times, and CV = 0.05 for all vehicle 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, whose values are LSL = 2 °C and USL = 7 °C in our analysis.Consequently, as indicated in Section 4.1, the index C pk has been adopted.In the sequel, the capability index C pk is computed for every route sequence following the Chen & Ding (2001) method described in Section 4.3, since it is computationally easy to implement and has shown to be superior when compared to other approximate PCI methods (Gonçalez & Werner, 2009).Although, in some manufacturing processes, values of C pk greater than 1.33 are adopted, a capable level, corresponding to a six-sigma formulation (see Table 1), seems more appropriate for this kind of process, since many external factors generate greater dispersions in the operational variables.Other authors, such Estrada-Flores & Eddy (2006), analysing a similar problem, have adopted the same C pk threshold.
The computer model (steps c and d, Section 2) was programmed in Free Pascal, version 2.6.4,to be further incorporated into a Delphi format.Applying the SA algorithm, the first stage involved 1,306 steps, and the second 10,815 steps, leading to the following optimal delivery sequence: Depot -8 -7 -4 -6 -9 -10 -11 -3 -12 -5 -2 -1 -Depot, with a total vehicle distance of 222.7 km, and with a satisfactory "capable" C pk value of 1.33.
The TSP alternative i.e., the route with minimum vehicle travel distance (204.1 km, in this application), yielded C pk = 0.35, well bellow the required standard level of 1.33.Taking the inverse TSP route, the result Figure 3 shows that, for this application and with a 95% confidence level, the four last visits indicated in Table 3, namely clients 12, 5, 2, and 1, should be assigned to another truck.In fact, a better approach to set the process within the specified temperature range should be the re-assignment of delivery stops per vehicle round, with the selection of smaller vehicles to perform the tasks defined with the aid of an optimization algorithm.
Although food temperature control largely influences microbial activities, allowing for an adequate surveillance of product quality along short and medium distance delivery routes, new risks of foodborne illnesses associated with micro-organisms are presently identified.One factor is the changing characteristics of the relevant micro-organisms.Second, the occurrence of changing production methodologies, in parallel with changes in the environment and an increase of the global trade of food stuffs (Havelaar et al., 2010).In parallel, it can be noticed in the literature substantial advances in the development of predictive microbiology, with efforts exploring a better balance between science and applications (McMeekin et al., 2008).Since temperature abuses in the distribution of refrigerated and frozen products are quite common in developing countries like Brazil (Pereira et al., 2010), its correct monitoring is naturally a first step in establishing a more comprehensive food perishability control, as is the case of the present work.Following this line further, the authors intend to research on the distribution of refrigerated products applying mathematical models describing the effects of environmental conditions on the growth and inactivation of microorganisms, in food distribution under more diverse logistics conditions.

Conclusions
The investigation presented in this paper has shown that the assignment of refrigerated vehicles to the distribution of perishable products is not a trivial process, requiring, in fact, a more detailed and integrated investigation of its thermal aspects.The use of process capability index evaluation over TTI data, together with a simulating annealing algorithm has produced an optimal constrained result in terms of a minimum-distance distribution route, but at same time respecting a pre-defined quality requirement, expressed through a C pk coefficient.travel times, including the line-haul segment and the local vehicle displacements, represented by the links from one served client to the next.For each stage, the simulation was replicated 1,000 times.The results are shown in Table 3 and Figure 3, where the mean C pk values are exhibited, as well as their lower and upper 95% significance levels.
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 fact has practical implications in real-life operations because different delays are common as, for example, unexpected 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 times, etc.And conversely, vehicle travelling times, both along the line-haul segments and along m 3 and m 4 being the second, third and fourth central moments of f (θ).
23) distance, but at same time keeping the TTI values within the required range.The TTI results from the SA model are represented by 394 temperature and time values, covering a typical tour.Figure 2 depicts the histogram of the TTI values obtained with the SA meta-heuristic.It can be seen that the temperature distribution is positively skewed (not normal), justifying the use of the C pk index (see Section 4).

Figure 2 .
Figure 2. Optimal SA route -histogram of internal vehicle temperature

Table 1 .
Quality status and C p values.

Table 2 .
Product quantity delivered per retailer client and mean unloading time.

Table 3 .
Optimal simulated annealing route: C pk values for each process stage.
Figure 3. C pk estimated values per routing stage