A Routing Model for the Distribution of Perishable Food in a Green Cold Chain

Abstract

In this research, we develop an extension of the stochastic routing model with a fixed capacity for the distribution of perishable products with a time window. We use theoretical probability distributions to model the life of transported products and travel times in the network. Our main objective is to maximize the probability of delivering products within the established deadline with a certain level of customer service. Our project is justified from the perspective of reducing the pollution caused by greenhouse gases generated in the process. To optimize the proposed model, we use a Generic Random Search Algorithm. Finally, we apply the idea to a real problem of designing strategies for the optimal management of perishable food distribution routes that involve a time window, the objective being to maximize the probability of meeting the time limit assigned to the route problem by reducing, in this way, the pollution generated by refrigerated transport.

1. Introduction

The need to transport products from an origin to a destination with one or more consumers requires the deployment of activities such as preparation, packaging, transportation, refrigeration, unloading, final storage and distribution to retailers. During the transportation of the product, the ideal storage, packaging and cooling conditions must be monitored to guarantee its quality. The process of avoiding the negative effects that the supply chain could have on the environment is called the green supply chain, which includes hazardous materials used during chain operations and incorporates environmental pollution, energy consumption and carbon emissions. A failed strategy in the design of the chain favors the decomposition of the product with the consequent loss of it and its equivalence in monetary units. In addition, there is pollution generated by poorly planned cargo transportation, poorly laid out routes, inadequate transportation and other factors that economically and environmentally affect the system associated with chain management.

The science of cooling consists of determining the conditions, physical, chemical, biological and engineering, under which low temperatures are applied in order to prepare the product for handling and travel and then return it to normal temperature without losing its organoleptic properties (odor, color, texture and flavor). Cryopreservation is the process of cooling and storing food, vegetables, fish, cells, tissues, organs, etc., at low temperatures or freezing them for future use. The need to transport cryopreserved products in a sustainable supply chain requires that it be efficient in aspects such as transportation techniques in refrigerated vehicles, freezing systems, stated routes and time windows. Currently, the efficiency of cryopreserved product chain management is remarkably improved through the use of emerging technologies such as Blockchain and Internet of Things (IoT) due to the impact on product greening. However, the factor with the greatest impact on the optimal management of the supply chain is the design of the route on time and at the minimum cost, considering the freshness of the transported products as the most important technological restriction of the plan. From this perspective, the most important issue is the design of the route that the transported shipload will have to follow to stay within the time window of the total freight.

In this research, we develop a stochastic mathematical programming model to obtain the optimal routing strategy based on the time window of the total freight since, regardless of emerging technologies, it is the lifetime of the products transported that defines the maximum duration of their transportation. Our project is based on defining the time window for loading a transportation system as a function of the minimum useful life associated with the products it transports. The novelty of our approach consists of integrating into a single real application model the characteristics of limited transport capacity, random time windows, random demands and routing times and the level of customer satisfaction included. In this regard, there are several authors who have addressed the problem with approaches similar to ours, among the most important of which are [1,2,3,4].

Our most important contributions to the problem are to consider, in a single routing model, the time window restrictions, the limitations on the carrying capacity of the vehicles, the randomness in demand and travel times, the clustering problem and the atomization of the products in primary and secondary routes. Above all, our objective function is expressed in terms of maximizing the probability that the travel time is less than the imposed time window. These types of functions are not common in the literature on the subject. Finally, we consider of utmost importance the development of a random search algorithm that satisfactorily approximates a solution to the model.

The organization of this document is as follows: in Section 2, we present a broad analysis of the economic importance of the cold chain and the consequences of its poor management in terms of the environmental pollution that this generates, especially in the form of greenhouse gases, which is used to justify the importance of the problem addressed; in Section 3, our development is formalized; in Section 4, we apply the model to a real problem of designing strategies for the optimal management of perishable food distribution routes, illustrating in detail the construction process and the results obtained; finally, in Section 5, we give the conclusions. A final appendix (Appendix A) accompanies this document.

2. Background

A supply chain is defined as the set of activities required for the delivery of goods or services to the consumer [5,6]. In a manufacturing line, the activities related to the supply chain define the processes necessary to convert raw materials or components into finished products or services. An essential element in the supply chain is logistics. This activity consists of the process of planning and executing transport, the efficient storage of goods and the correct distribution of the products between customers and suppliers [7,8]. In this case, transport is the central problem for the design of the chain; this is where routes, loads, transport typology, operators, travel times, transport costs and time windows have to be defined. The importance of each of these issues is as described below:

• Routes. This considers the activity consisting of choosing the paths that the product distribution units must travel, considering distances, types of highways, their physical condition, highways, tolls, dirt roads, diversions, traffic regulations, assistance centers, routes for dairy (the milk road) and primary and secondary routes.
• Volumetry. This includes the calculation of areas and volumes available inside the transports, also considering the wasted spaces and the logistics of product accommodation inside it. Normally, the Last Input First Output (LIFO) discipline is used for the access and exit of containers in each transport, although this depends on the design of the vehicle, the container and the needs of the transport company, as well as the requirements of the clients.
• Type of transport. This concept includes land, sea, air, intermodal, river transport, pipeline transport and Ro-Ro transport (the acronym for “Roll on-Roll off” that refers to the system by which a ship transports cargo on wheels, mainly cars or trucks). In relation to land transport, it must be specified if the vehicle is a dry box, refrigerated dry box, wet cargo transport, platform, curtain, simple container, refrigerated container, heat truck or tank truck. In other types of transportation such as sea and air, there are other technological restrictions that must be addressed when shipping a load.
• Operators. This refers to the designation of the best transport driver based on his expertise, experience, knowledge of the route, travel times, rest time, time to eat and time to wash. The problem consists of assigning an operator to a route, a transport or a load.
• Travel times. This includes from the moment the vehicle forms in a port, bay or platform to receive the load until the moment its final point of travel arrives and the total load has been unloaded. This includes areas of spraying or primary routes and the handling of secondary routes. The geographical place where the full load is divided into smaller loads to be distributed to the retail centers is called the shoveling zone. The route followed for this activity, using smaller vehicles, is usually called a secondary route. Similarly, a primary route is defined as the route that a transport with a large load capacity (usually more than one ton) must cover to take the product from the distribution center (DC) to the shoveling nodes. In both cases, not only is travel time important, but also the rest time, cleanliness of the vehicle operator and his food.
• Time window. This defines the interval that exists from the moment the transport load is made until the moment the product is delivered to the final customer; that is, the last link of the chain. This includes the times mentioned in the previous item. The time window helps determine the minimum time a carrier must take to deliver the product before it degrades. The product that must be transported with refrigeration in order to extend its lifetime deserves special mention.

2.1. Product Spoilage

A perishable product (PP) is defined as one that, by natural consequence, deteriorates in terms of its quality and/or organoleptic properties. This may be due to its biological and/or physical–chemical properties, for example, in the absence of adequate storage conditions. For this reason, packaging, loading, transportation and unloading, and where appropriate, refrigeration, must be performed under controlled conditions in order to minimize loss or damage. Examples of these products are fish, fruit, vegetables, potatoes, plants, bread, meat (of all types), hunted animals, butter, eggs, milk, cheese, birds, small animals, etc. Storage for these generally requires temperature control and proper treatment for transportation [9].

The concept of deterioration refers to the process of diminishing the quality, function and/or condition of the product in relation to its original attributes. As an example, damage to fruits and vegetables can occur in the form of rot, putrefaction, expiration, etc. In this case, the damage is due to the product continuing to metabolize and respire after harvest. Its loss is based on the fact that oxygen from the atmosphere is absorbed, producing heat, carbon dioxide and even ethylene gas. Here, the respiration rate is defined as the amount of $O_2$ consumed or $C{O}_2$ produced by the product, per mass unit and per time unit. In this case, the respiration rate has been identified as one of the indicators of physiological stress. Therefore, to keep products fresh and prolong their shelf life, the respiration rate must be reduced without affecting their organoleptic properties [10].

On the other hand, there are two factors of great importance to consider in the PP supply chain. The first one is the economic factor associated with the harvest, storage and transport of the product from its origin to the final consumer. Secondly, there is the environmental impact derived from the management of the supply chain associated with its distribution process. More generally, the abuse of temperature in the Food Cold Chain (FCC) is a problem of crucial interest for the production and the logistics [11], as well as the systematization of the management of the delivery process of perishable cargo and the development of new methods to improve its efficiency.

The transport of PPs requires greater attention to maintain the quality of the cargo. This is achieved through careful preparation and strict adherence to numerous requirements for loading, location, packaging, delivery time, etc. The foregoing requires high demands on the transport organization process in terms of the coherence and coordination of the actions of the actors involved. In addition, it requires compliance with standards and quality planning of both the process and the delivery [12]. Finally, the strict regulation of quality in the supply chain and its logistics process, as well as its economic indicators, are given in the SCOR model (Supply Chain Operations Reference) [13]. In this reference, the main standards that producers, carriers, warehouses, sales and distribution centers must comply with are defined.

2.1.1. The Environmental Impact of PPs

The useful life, overproduction, storage and especially the transport of PPs are factors that have a strong effect on the environment. Therefore, there is a requirement to integrate sustainability with the design of the supply network. PPs generally have a short useful life, and their storage must be constantly controlled to prevent expired waste and its impact on the environment [14]. The importance of good design in the supply chain affects economic action and environmental pollution. Some aspects related to pollution are discussed in detail below.

The design of the PP supply chain must consider its expiration time. Therefore, the time windows play an important role in the topology of the product distribution network. Some of the most relevant factors associated with the environmental contamination generated by the use of transport and the handling and cryopreservation of PPs are the following:

• Emission by polluting particles from the combustion of transportation engines and greenhouse gas (GGE) emissions.
• Contamination generated by the refrigeration systems in the plant and during handling.
• Contamination generated by packaging.
• Pollution generated by the expired product.

In general, environmental impact is defined as the direct effect of socioeconomic activities and the natural consequences on the components of the environment. It can be said that, for the most part, it is due to the activities of the human being. In [15], it is recognized that this phenomenon has three components:

• The environmental stressor: the pollution and noise that can be measured in terms of tons of transported product divided by tons of pollutants emitted into the atmosphere.
• The spatial pattern of the distribution of transported goods: The mode of transport and the total amount of stress placed on the environment. This depends on the volume of the products transported and the distance traveled.
• The environmental impact: the nature of the environment; for example, the characteristics of the physical ecosystem and human density.

Regarding the previous list, the air pollution generated by transport is the main cause of its poor quality and directly affects the health of living beings. Pollutants that cause poor air quality include particulate matter, nitrogen oxides ($N{O}_x$) and volatile organic compounds (VOC). In general, the transport sector is responsible for the emission of approximately the following figures of pollutants [16]:

• More than 55% of $N{O}_x$ emissions.
• Less than 10% of VOC emissions.
• Less than 10% of emissions of very small particles (organic chemicals, dust, soot and metals) suspended in the air that have a diameter of less than 2.5 microns. Also, there are small solid or liquid particles of dust, ashes, soot, metal, cement or pollen, whose aerodynamic diameter is less than 10 micrometers.

The following percentages show the distribution of greenhouse gas emissions worldwide during the year 2018 [17]:

• Electric power generation: 32%.
• Transport: 17%.
• Manufacturing and construction industry: 13%.
• Agriculture: 12%.
• Industrial processes: 5.9%.
• Fugitive emissions (spurious leaks in industrial areas and/or clusters of companies such as gases used in refrigeration systems and others that come from equipment such as valves, pumps, pipes, etc.): 5.9%.
• Residential areas: 5.9%.
• Waste of all kinds: 3.3%.
• Other combustion sources: 3%.
• Land use, land use change and forestry: 2.8%.

The above amounts show only part of the environmental pollution problem generated by the emission of gases and particles from transportation. As a consequence, pollution in the US alone is to blame for approximately 40,000 premature deaths, 34,000 hospitalizations and 4.8 million lost work days. Therefore, the US Environmental Protection Agency (EPA) expects that emissions of hazardous air pollutants from mobile sources will be reduced by 80% by the year 2030.

Other important factors associated with pollution are the following:

• Operational oil contamination.
• Solid waste disposal.
• Accidental spills.
• The construction and maintenance of ports and canals.
• Pollutants from the aeronautical and railway industries, due to the use of ducts and pipes.

With respect to the contamination generated by the refrigeration systems in the plant and during handling, the effect of air and water pollution on industrial refrigeration systems has been well documented. Air contamination reduces the efficiency of the refrigeration system and increases electrical costs due to the increased compressor discharge pressure caused by the presence of air.

The problem with the use of refrigerants is that they are known to have a negative effect on the environment. The above is due to its contribution to global warming and the depletion of the ozone layer by what are called greenhouse gases (GHGs). These, like carbon dioxide and emissions from other types of refrigerants, contribute to global warming by absorbing infrared radiation and retaining it in the atmosphere. For example, Chlorofluorocarbons (CFC), Hydrochlorofluorocarbons (HCFC) and HFC Hydrochlorofluorocarbons (compounds made up of hydrogen, fluorine and carbon).

The impact of air conditioning and refrigeration systems on stratospheric ozone is mainly related to the release of refrigerants that deplete it. This is one of the main causes of global warming since it contributes through the release of refrigerants and GHGs. In this case, the electrical energy necessary to operate the refrigeration systems also has a significantly greater warming impact. Hence, the importance of phasing out Hydro-Fluoro-Carbon (HFC)-based refrigerants with less efficient options will increase net GHG emissions; the same conclusion applies to Per-Fluoro-Carbons (PFC), although they are used less frequently as refrigerants [18].

The increase in atmospheric concentrations of CFC has accounted for about 24% of the direct increase in the radiative forcing of greenhouse gases during the last decade. A decrease in the amount of stratospheric ozone has been observed and is believed to be related to the increase in stratospheric chlorine from CFCs (and, to a lesser extent, from other man-made compounds containing chlorine and bromine) [19].

The food distribution process using refrigerated systems is responsible for a significant loss of quality in perishable products and a huge source of environmental pollution in addition to waste from cooling systems (air conditioners, freezers, refrigerators, chillers and dehumidifiers) [20]. Due to the above, one of the main lines of current research in the area of Engineering in Cold Technology (ECT) consists of the search for new ideal refrigerants that are less harmful to ozone; that is, compounds that are free of these problems but that, in turn, provide safety, stability, compatibility, cost and similar burdens. Some of these trends incorporate $C{O}_2$ as new refrigerant possibilities associated with thermoelectric and ejector-compression cooling system technologies; see, for example, [21].

Waste from cooling systems (air conditioners, freezers, refrigerators, chillers, dehumidifiers, etc.) is another major source of air pollution. In the case of a refrigerator or freezer, they generate an environmental impact made up of the following compounds: insulating foam (9 kg), refrigerant gas (125 g), contaminated oil (250 g), mercury (above 2 g), plastics, metals and glass (76.5 kg) and carbon dioxide (above 2.5 tons) [22]. In relation to the operation in processing plants, aerosols are also common sources of pollution after industrial processes such as the pasteurization of dairy products. Little is known about the degree to which biological aerosols contaminate pasteurized products and the environment in general; however, evidence indicates that the air within a packaging area is a critical control point for both pathogens and spoilage microorganisms.

Regarding the manufacturing of dairy products, the microorganisms involved in the product are often damaged due to the stresses of the aerosol state and consequently may not grow on selective media; for example, in cheese production, yogurt, etc. Aerosols are generated within the dairy by worker activity, sink and floor drains, water spray and air conditioning systems [23]. Similar situations are experienced in processing plants for meat products and their derivatives, the bread industry, sausages, etc. The effect caused by the oil used in the systems of refrigeration relates to its use as hydraulic control, functional fluid and lubricating oil in refrigeration compressors under the influence of a refrigerant. Cooling systems are of two types: synthetic oil and mineral oil. The first one has more duration, and the second one is used in industrial applications such as air conditioning units for commercial buildings [24].

Relating to the contamination generated by packaging, this is the wrapping or bottling of the products to protect them from being damaged during handling, transport and storage. Its function is to keep them safe and marketable by also helping to identify, describe and promote them. The most common materials used in product packaging are the following:

• Rigid plastic containers or PET (polyethylene terephthalate) or high-density polyethylene plastic (HDDE).
• Paper.
• Cardboard.
• Cardboard/Fiberboard.
• Aluminum.
• Glass.
• Expandable polystyrene (styrofoam).
• Flexible plastic containers.

Polyethylene terephthalate (PET) is a polymer not classified as a dangerous substance according to Regulation (EC) No 1272/2008 (CLP). PET is not classified as persistent, bioaccumulative or toxic (PBT). However, its use generates serious problems for its disposal since its recycling requires large amounts of energy with conventional sources, which depletes natural resources and generates environmental degradation [25]. PET plastic is the most widely used in the production of single-use plastic-based water containers. It does not contain BPA, but PET is also associated with many of the same health risks, including delayed growth, reproductive problems, low energy levels, body balance problems and an inability to process stress [26]. An interesting figure is that single-use plastic bags take approximately two decades to degrade. In contrast, plastic water bottles made with polyethylene terephthalate are estimated to take approximately 450 years to fully decompose.

For its part, paper alone represents approximately 40% of all waste around the world. That adds up to around 71.6 million tons per year in the US alone. Paper waste is a huge problem in terms of the devastating impact on the world due to deforestation. Paper production uses up to 40% of all the world’s wood. The process of paper manufacturing releases nitrogen dioxide, sulfur dioxide and carbon dioxide into the air, contributing to pollution such as acid rain and greenhouse gases [27].

Cardboard is a heavy, stiff paper used to make boxes used for packaging. It is made of several layers of thick paper to increase its rigidity and strength in order to protect the items to be stored. Cardboard drums are even used to transport dangerous chemicals, products pharmaceuticals and hazardous waste. Cardboard is biodegradable, produces methane (a GHG) as it breaks down and typically ends up in a landfill, increasing the amount of methane released into the atmosphere.

The contamination produced by aluminum occurs mainly during its extraction process. Some toxic pollutants produced during this are dioxins, furans, hexafluoroethane, tetrafluoromethane, fluoride (gases and particles), hydrocarbons polycyclic aromatics (PAHs), mercury and benzopyrene.

Expanded polystyrene (EPS) is a white foam plastic material produced from solid beads of polystyrene. It is mainly used for packaging, insulation, etc. It is a closed-cell rigid foam material produced from (a) styrene—which forms the cellular structure and (b) pentane, which is used as a blowing agent. Plastic can be degraded into microparticles (MP) < 5000 nanometers in diameter and then into nanoparticles (NP) < 100 nanometers in diameter. NP have been detected in air, soil, water and sludge. The use and handling of this substance entails health risks; for example, polystyrene nanoparticles can penetrate organisms, accumulate throughout the food chain and are also surrounded by a crown of proteins that allows them to penetrate membranes and are highly toxic.

Depending on the cell type, NP can be transported via pinocytosis, phagocytosis or passively transported. Currently, there are no studies indicating the carcinogenic potential of NP. On the other hand, the PS (styrene) monomer was classified by the International Agency for Research on Cancer (IARC) as a potentially carcinogenic substance (carcinogenicity class B2) [28].

Finally, the contamination generated by the expired product is considered from cosmetics to food. In the second case, food and ingredient manufacturers often include production dates, best-before dates and/or spoilage information. Foodborne contamination is generally defined as food that spoils or is contaminated because it contains microorganisms, such as bacteria or parasites, or toxic substances that make it unsuitable for consumption. A food contaminant can be biological, chemical or physical in nature, the former being the most common. The four main types of contamination considered in this document are chemical, microbial, physical and allergenic. All foods are at risk of contamination with these four types. Even canned foods expire; for example, canned foods with a high acid content, such as tomatoes and fruits, can expire after 12 to 18 months. Low-acid canned foods, such as meats and vegetables, can stay fresh for 2 to 5 years. However, if the containers are rusty, dented, swollen or otherwise damaged, they indicate that the food is unsafe or may be contaminated. In relation to frozen foods, they do not expire.

When food expires, the vast majority of it is thrown away. This is not just an obvious problem for food safety: it is a huge environmental problem. In addition to the amount of land and water required to produce all the food that is never used, the Food and Agriculture Organization of the United Nations has estimated that the carbon footprint associated with food waste worldwide each year is more than 3 billion tons of carbon dioxide equivalent. The FAO reported that around 1.6 billion tonnes of food is wasted globally each year, thanks to a variety of reasons at points along the supply chain. In the US, research has suggested that up to 40% of the nation’s food supply could end up wasting away [29]. Other PPs that deserve special attention due to their environmental impact are cigarettes, food wrappers/containers, beverage bottles, plastic bags, caps/lids, cups, plates, forks, knives, spoons, straws/stirrers, glass beverage bottles, beverage cans, paper bags, etc. [30,31].

2.1.2. The Design of the Supply Chain in Perishable Products: A Short Bibliography Review

The design of supply chains for products of all kinds is a classic problem of Industrial and Manufacturing Engineering as well as Cold Technology Engineering in order to provide the producer and customer with the optimal means for the manufacturing, packaging, loading, transport and download of the product with time windows. The concept of optimality refers to the delivery time of the product (lead time) in the quantity and quality requested and fully complying with the regulations established by the SCOR model, which is the reference framework for supply-chain operations. Design problems also involve the macro- and microlocation of manufacturing centers, storage, spraying areas, loading and unloading areas and delivery criteria in cities, ports and other terminals where the product must be stored. This includes some of the characteristics of transport, routes and packaging systems involving, in addition, the refrigeration requirements. Practical applications of these models can be seen in [32,33].

The literature in this regard is abundant in terms of information related to models for the design and operation of supply chains. Little attention has been paid to the models associated with perishable food and even less attention has been paid to products that require cryopreservation. Regarding production planning and the logistics of perishable products, the investigation is scarcer, in particular under the approach in which this point of view is related to time windows, which will be addressed. For example, in [34], an integrated location-inventory-routing model for perishable products in an emerging market is developed. The analysis includes the cost, freshness and carbon-emission factors. The authors propose a multitarget planning model and set constraints based on real location-inventory-routing situations.

In [35], approximation methods are applied to estimate and compare the total logistics cost of supply-network designs under various business conditions, such as variations in demand, changes in costs and changes in production policies. The model is applied in a real case and evaluates the trade-offs between five different network designs for the supply of highly perishable foods from a single regional supplier. Similarly, in [36] a multi-objective mathematical programming model is developed to optimize the cost, energy consumption and traffic congestion associated with supply chain operations. In [37], the authors elaborate a systematic literature review of articles focused on sustainable supply-chain management in global supply chains.

In [38], a new idea is developed to study the use of environmental indicators to measure profitability. The objective of the authors was to understand the impact of a group of categories (on the installation, location and capacity of the installation, selection of suppliers, technology used, definition of the transport network, supply planning and product recovery) in the decisions of the chain. The model was applied to a real case study.

In the case of systems with refrigeration, in [36], the authors develop a multiobjective mathematical programming model to optimize the cost, energy consumption and traffic congestion associated with the operations of a supply chain. The important contribution of the model is very close to the objective of our development since we analyze the useful life of the product through a Weibull-type random variable assuming that the expiration date of the food is affected by the use of the vehicle’s refrigerator, which is considered a decision variable. The model includes different classes of vehicles and multiple types of products. A dairy-supply-chain case is investigated and sustainability interrelationships are studied. A similar study is the work presented in [39]. Here, the degradation process of perishable foods is studied and the optimal temperature of the food chain is determined, as well as the optimal price to maximize the benefit of the channel. The optimum temperature is correlated with the profit proposed for its sale.

An interesting investigation is found in [40]. This article aims to determine the principles of the sustainable distribution of perishable goods and examine the current state and plans for the application of its principles in business practice on the example of port cold storage. The authors propose market research to identify the main directions of activities carried out by port cold stores in the field of the sustainable distribution of perishable foods.

In the field of mathematical modeling similar to the one developed here, there are important contributions in this regard; for example, in [41], a methodology is proposed for the rapid loading of products in an electric vehicle assisted by a service vehicle (drone). In [42], the authors address the stochastic dynamic vehicle routing (SDVR) problem inherent to urban logistics and from a theoretical perspective, and ref. [43] addresses the problem of the minimum spanning tree, a fundamental tool in the analysis and construction of clusters in logistics engineering.

A very important contribution to the topic is found in [44]: a vehicle-routing problem with time windows (VRPTW) with compatibility matching constraints and the total completion time as the objective function is proposed, with applications in home healthcare routing and scheduling. An approach close to ours is found in [45]. Here, the problem of recycling waste products through closed-loop logistics networks is analyzed by using a fuzzy programming model to minimize the total network cost and the sum of carbon rewards and penalties when selecting recycling locations, facilities and transportation routes among network nodes. Other representative models of the cold-supply-chain-management problem are given in [46,47,48].

Despite the technological advances applied to solving stochastic models, there are still several areas of opportunity to be addressed in future research. For a comprehensive review of the methods and models associated with the routing model, see [49].

3. Focusing on the Routing Problem with a Time Window: Setting Up the Mathematical Model

Once the importance of our offering has been justified, we will now focus on developing the mathematical model of optimization. The classic model used in the design of distribution networks is called the routing model, a combinatorial optimization problem. Depending on the characteristics of the topology of the distribution network, it is common to design some instance when deterministic information is available. The main components of the model (regardless of their topology) are:

• A product distribution center.
• A set of customers who demand a quantity of products with a certain demand (not necessarily known).
• A fleet of available vehicles with a known transportation capacity.

Our assumptions in this research consider the existence of primary and secondary routes as well as spray nodes along the trajectories inherent to the distribution network. We will also assume that the shelf life of a perishable product can be uniquely modeled through a theoretical or empirical probability distribution function, and the same case applies to the estimation of demand and travel times along the network. The problem is to plan the delivery of the requested products at the minimum cost or time considering that vehicles leave and return to the same warehouse with certain operational restrictions. In this, we only have to deal with the computational complexity of such an instance due to the number of variables involved in the model.

With respect to the above, the most common variants of the vehicle-routing models are explained below [50,51,52,53,54,55,56]:

• Vehicle routing with limited capacity.
• Capacitate symmetrical.
• Vehicle routing problem (VRP) with time windows.
• VRP with transshipment warehouses.
• VRP with stochastic demand.
• VRP with regular deliveries.
• Trained VRP with delivery and pickup.
• VRP with shoveling centers.
• VRP with open routes.

Similarly, the topology of the networks used in these models and/or their variants have the same characteristics; for example, a central warehouse with two or more distribution routes (Figure 1) or a central warehouse with several shoveling centers (also called spraying) from where the products are distributed to peripheral clusters (Figure 2).

In this project, we are interested in complex networks with intermediate stores and a double routing system. The first routing component is called the primary route, and its purpose is to decide the path from the center deposit to the intermediate warehouses and from there to the respective stems that communicate with the system cluster. The secondary route must define the route that will be made with light transport from the shoveling node to the consumption centers that make up the cluster. The working hypothesis of this document is that the topology of the distribution network is known and fixed; that is, the clustering process and the sequencing of the routes (the traveling agent problem) to be taken are ruled out in this model. The research question is: can the delivery of a perishable product be guaranteed with time window limitations and transportation and loading information in a complex vehicle-routing model while meeting a level of customer service? The solution to this question will be approached from the perspective of maximizing the probability of the delivery of merchandise within the imposed sales time.

3.1. Materials and Methods

For the construction of this model, a stochastic mathematical programming model is developed in order to select the minimum route of the system’s transport network. The model has the following structure [57]:

$$\begin{array}{c}{\mathrm{Minimize}}_x{g}_0\left(x\right)={\mathrm{Minimize}}_x\int {\psi }_0(x,y)f_Y(y,x)dy={\mathrm{Minimize}}_x\mathbf{E}\left[{\psi }_0(x,Y)\right]\hfill \\ \mathrm{Subject} \mathrm{to}\hfill \\ g_j\left(x\right)=\int {\psi }_j(x,y)f_Y(y,x)dy=E\left[{\psi }_j(x,Y)\right]\le 0,j=1,\dots ,M,x\in X\subset {\mathbb{R}}^n\hfill \end{array}$$

(1)

where x and y are n and m dimensional vectors, respectively, and $f_Y(y,x)$ is a probability density function (usually unknown) of a random variable Y that depends on a vector of parameters x. A great variety of applications and optimization methods exist in the literature for different variants of this problem.

Model Construction

The construction of this model is based on the topology of Figure 2. The first assumption in this document is that the lifetimes $T_j$ of transported products can be modeled through a probability distribution function (theoretical or empirical). For the above, the lifetime of the product i transported can be viewed as a set of random variables. Therefore, the strategy used consists of building a network similar to the one shown in Figure 2 and developing a stochastic mathematical programming model associated with it. The second part of this writing consists of showing an instance associated with the model and the corresponding optimization method used. The construction of the proposed model will be performed step by step to illustrate its nature. These ideas are developed below.

Let $T_1,T_2,...,T_{\theta }$ be a set of independent random variables defined over a probability space $(\Omega ,\mathcal{F},\mathcal{S})$ (not necessarily with the same distribution) that denote the lifetime of the $\theta$ transported products. Then, the required time window of the total transport load $W_t$ is given by

$$W_t=min\{T_1,T_2,...,T_{\theta }\}$$

(2)

Analogously, we will also assume that, for every $T_i$, there is a function $F_{T_i}\left(t\right)$ such that

$$F_{T_i}\left(t\right)={\int }_0^t{f}_{T_i}\left(t^{\prime }\right)d{t}^{\prime },t\ge 0$$

(3)

where

$$\mathbf{E}\left[T_i\right]={\int }_0^{\infty }td{F}_{T_i}\left(t\right)<\infty ,\mathrm{and}\mathbf{Var}\left[T\right]=\mathbf{E}\left[T^2\right]-{\left(\mathbf{E}\left[T\right]\right)}^2<\infty ,$$

(4)

where $\mathbf{E}$ is the expectation operator. By the continuity hypothesis of T, $F_{T_i}\left(t\right)$ is a nondecreasing function, and therefore, there is the inverse function ${\xi }_i=F_{T_i}^{-1}\left(u\right)$ defined for any value of u in $[0,1]$ such that

$${\xi }_i=F_{T_i}^{-1}\left(u\right)=inf\{t:F_T\left(t\right)\ge u\},0\le u\le 1$$

(5)

According to Figure 2, the vehicular routing must cover two well-differentiated routes, which we will call the primary route and the secondary route. The primary route is intended to route the goods from the producing area to the shoveling nodes, moving them by means of large-tonnage transport. For each shoveling node, there is an associated cluster where smaller vehicles distribute the merchandise in smaller quantities to final customers. This route is called secondary.

In the first part of this development, we will address the travel times in the network. Assume that the primary route is a set S consisting of K shoveling nodes. In addition, each shoveling node $k\in K$ has an associated set $N_k$ (cluster) of shoveling nodes whose cardinality is $\phi$. Formally,

$$S=\{s_1,s_2,\dots ,s_K\}$$

(6)

$$N_k=\{n_1,n_2,\dots ,n_{\phi }\},k=1,2,\dots ,K$$

(7)

Let ${\mathcal{T}}_0$ and ${\mathcal{T}}_2$ be the travel times from the production center to the packing warehouse and from the packing warehouse to the first shoveling node, respectively. Also, let ${\mathcal{T}}_1$ be the occupation time of the product in the packing warehouse. Let us also define the variables ${\mathcal{T}}_{jk}$ as the travel time from shoveling node j to shoveling node k, $k>j$, $k=2,\dots ,K$. This means that the total travel time of the primary route is given by

$${\mathcal{T}}_f={\mathcal{T}}_1+\sum _{(s_j<s_k)\in S}\left({\mathcal{T}}_{jk}+{\mathfrak{d}}_{jk}\right),j=1,\dots ,K,k=2,\dots ,K+1$$

(8)

where ${\mathcal{T}}_1={\mathcal{T}}_0+{\mathcal{T}}_1+{\mathcal{T}}_2$ is the total time required to reach the first shovel node and ${\mathfrak{d}}_k$ is a random variable that denotes the loading time + the unloading time of the products in and from transportation in the shovel node k.

Similarly, let ${\tau }_{ij}$ be the random travel time between node i and node j with $n_i$ and $n_j$∈$N_k$; then, the total travel time ${\tau }_k$ on a secondary route k is given by

$${\tau }_k=\sum _{(n_i<n_j)\in N_k}{\tau }_{ij},k=2,\dots ,K$$

(9)

From the above, the following facts can be deduced

• The total time required to reach the first shovel node ($s_1$) is

  $${\mathcal{T}}_1={\mathcal{T}}_0+{\mathcal{T}}_1+{\mathcal{T}}_2$$

  (10)
• The total time required to complete the traversal of the secondary path associated with the first shovel node is

  $${\omega }_1={\mathcal{T}}_1+{\tau }_1$$

  (11)
• The total time required to reach the second shovel node $\left(s_2\right)$ is

  $${\mathcal{T}}_2={\mathcal{T}}_1+{\mathcal{T}}_{12}+{\mathfrak{d}}_1$$

  (12)
• The total time required to complete the traversal of the secondary path associated with the first shovel node

  $${\omega }_2={\mathcal{T}}_2+{\tau }_2$$

  (13)
• It is satisfied, in general, that

  $${\mathcal{T}}_K={\mathcal{T}}_{K-1}+({\mathcal{T}}_{K-1,K}+{\mathfrak{d}}_{K-1})$$

  (14)

  $${\omega }_K={\mathcal{T}}_K+{\tau }_K$$

  (15)

  represent the total travel time of the main route and the total travel time of the k-th cluster.

Regarding the storage capacity of the packing warehouse ($Ca{p}_m$), we will assume that it is finite. Similarly, the existence of m cargo trucks will be assumed to take the product to the shoveling nodes. Then, to move $\theta$ products in any of the M transports, it is satisfied that

$$\sum _{r=1}^{\theta }{v}_r{x}_r\le Ca{p}_m,m=1,2,\dots ,M$$

(16)

where $v_r$ is the volume unit (in cubic meters per piece) of the r-th product $x_r$ (measured in pieces) and $Ca{p}_m$ is the capacity (in cubic meters) of the m-th vehicle.

In addition, all the products transported in each period must satisfy the consolidated demand D (in pieces) of all the clusters that are served:

$$\sum _{m=1}^{M}Ca{p}_m\ge \sum _{s\in S}{D}_s=\sum _{s\in S}\sum _{n\in N_k}{d}_{sn}=\mathcal{D}$$

(17)

Furthermore, in order to meet a customer service level $(1-\alpha )$ with $\alpha \in [0.01,0.05]$, we have that

$$\mathbf{P}\left(\sum _{m=1}^{M}Ca{p}_m\ge \sum _{s\in S}\sum _{n\in N_k}{d}_{sn}\right)=1-\alpha ,$$

(18)

Then, for each time a shipment is made, the following stochastic mathematical programming problem must be solved:

$$\mathrm{Maximize}Z=\mathbf{P}(W_t\le t^{\prime })$$

(19)

subject to:

$${\mathcal{T}}_K={\mathcal{T}}_{K-1}+({\mathcal{T}}_{K-1,K}+{\mathfrak{d}}_{K-1})$$

(20)

$${\mathcal{T}}_K+{\tau }_K\le t^{\prime }$$

(21)

$$\sum _{s\in S}\sum _{n\in N_k}{d}_{sn}=\mathcal{D}$$

(22)

$$\mathbf{P}\left(\sum _{m=1}^{M}Ca{p}_s\ge \mathcal{D}\right)=1-\alpha ,$$

(23)

$$\mathcal{D}\sim \mathcal{N}({\mu }_{\mathcal{D}},{\sigma }_{\mathcal{D}}^2),T_i\sim f_{T_i}\left(t_i\right),t^{\prime }>0.$$

(24)

3.2. Solution Strategy

The strategy followed to solve the proposed model consists of finding an equivalent model that is simpler to solve. Regarding the objective function, notice that for $t^{\prime }<\infty$, we have

$$\begin{array}{ccc}\mathbf{P}(W_t\le t^{\prime })=P\left[\min \{T_1,T_2,\dots ,T_{\theta }\}\le t^{\prime }\right]& =\hfill & 1-P(T_1>t^{\prime },\dots ,T_{\theta }>t^{\prime })=\hfill \\ 1-\left[[1-F_{T_1}\left(t^{\prime }\right)][1-F_{T_2}\left(t^{\prime }\right)]\cdots [1-F_{T_{\theta }}\left(t^{\prime }\right)]\right]& =\hfill & 1-\left[{\Pi }_{i=1}^{\theta }[1-F_{T_i}\left(t^{\prime }\right)]\right].\hfill \end{array}$$

(25)

In particular, when the random variables are identically distributed, the above expression simplifies to

$$\mathbf{P}(W_t\le t^{\prime })=1-{\left[1-F_T\left(t^{\prime }\right)\right]}^{\theta }.$$

(26)

Regarding the level of customer service, Equation (18) is equivalent to

$$\mathbf{P}\left(\Xi \le \varrho \right)={\int }_0^{\varrho }{h}_{\Xi }\left(\xi \right)d\xi =1-\alpha ,$$

(27)

where $\Xi$ is the random variable that represents the consolidated demand, which, for simplicity, we will assume to be normally distributed; that is,

$$\Xi \sim \mathcal{N}(\mu ,{\sigma }^2),\mathrm{with}\widehat{\mu }\approx \sum _{s\in S}\sum _{n\in N_k}{d}_{sn},$$

(28)

Using the following equivalence in density functions,

$$\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{\Xi -{\mu }_{\Xi }}{{\sigma }_{\Xi }}\right)}^2}\equiv \frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{z}^2},z=\frac{\Xi -{\mu }_{\Xi }}{{\sigma }_{\Xi }}$$

we have

$$\mathbf{P}\left(\Xi \le \varrho \right)={\int }_{-\infty }^{\varrho }\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{1}{2}{\left(\frac{\Xi -{\mu }_{\Xi }}{{\sigma }_{\Xi }}\right)}^2}d\Xi \equiv \Phi \left(\varrho \right)=1-\alpha ,$$

(29)

where $\Phi \left(\varrho \right)$ is the standard normal distribution $\mathcal{N}(0,1)$ evaluated at the point

$$\varrho =\sum _{s\in S}\sum _{n\in N_k}{d}_{sn}$$

Therefore, Equation (29) can be solved by using

$$\Xi ={\Phi }^{-1}\left(\varrho \right)={\sigma }_{\Xi }z+{\mu }_{\Xi }$$

(30)

or by using the error function as follows

$${\Phi }^{-1}\left(\varrho \right)=(2\varrho -1)\left(er{f}^{-1}\left(\varrho \right)\right)\sqrt{2},erf\left(\varrho \right)=\frac{1}{\sqrt{2\pi }}{\int }_0^{\varrho }{e}^{-u}du$$

(31)

3.3. Density Functions Used in the Model

The lifetime of a PP is never constant since variables intervene, such as the plant-processing time, packaging time, initial quality of the product, process used for cooling and time required to distribute it. Therefore, product lifetimes will be evaluated based on random variables with known densities. Thus, given the nature of our approach, we only used stochastic models in terms of known density functions by using proven methods to create the lotteries inherent in the process corresponding to the required samples. For example, the use of this approach to model the lifetime of human beings is kindly illustrated in [58]. In particular, we used a family of exponential distributions to model the lifetime of the products studied [59,60]. To represent the life of foods, fruits and vegetables, we use similar ideas similar to those shown in [61,62,63].

The literature on this matter is abundant and provides us with this approach, which is quite appropriate for realistic modeling. Due to their nutritional importance and physical appearance, the main quality factors (appearance, texture and flavor) are the organoleptic and sensory properties since they are perceived by the senses. Some of the variables measured during the process were [64]:

• The size, shape and integrity (refers to the degree of the whole and broken pieces).
• Color and shine.
• Consistency as an attribute of textural quality; food texture can be reduced to measurements of resistance to force. It can also be detected by sensory means such as fingers, the tongue, the palate or teeth.
• Texture changes such as softness, hardness, crystallization and more.
• Flavor evaluated subjectively through sensory sampling.
• Flavor panels. Some flavor-providing substances can be measured chemically or physically with other instruments. Some examples are salt, sugar and acid.
• Nutritional quality, sanitary quality and storage quality: Nutritional quality can frequently be assessed through chemical or instrumental analyses of specific nutrients. Sanitary quality was generally obtained through measurement and counts of bacteria, yeast, mold and insect fragments, as well as sediment levels. In relation to storage quality, this property is measured through the maintenance of quality, or storage stability is measured in storage and handling conditions that are configured to simulate or exceed to some extent the conditions that the product is expected to encounter in its normal distribution and use.

For the above, a variable sample of five typical perishable products (fresh chicken; fresh fish; raw sausages, turkey and pork; fresh vegetables; and eggs) is considered, which were subjected to quality tests until any of the failures inherent to the product appeared ($t^{\prime }$).

Table 1. Characteristics and sample sizes $\left(Ss\right)$ used to build the instance.

Product	Packaging	Packaging	Presentation	${\mathit{t}}^{\prime }$ (in h)	$\mathit{S}\mathit{s}$
	Conditions	Temperature	in kg
Fresh chicken	fillets	$-18°\mathrm{C}$	0.250 each	17.82	1000
Fresh fish	fillets	$-18°\mathrm{C}$	0.250 each	46.58	1000
Raw sausages	bottles	$-19°\mathrm{C}$	0.445 each	65.80	1000
Turkey and pork	fillets	$-18°\mathrm{C}$	0.250 each	62.43	1000
Fresh vegetables	package	$+14°\mathrm{C}$	0.125 each	32.08	1500
Eggs	piece	$+10°\mathrm{C}$	0.125 each	34.12	1500

shows the characteristics and sample sizes $\left(Ss\right)$ used in the project.

Once the sample values were obtained, the data were adjusted to some probability distribution. In this case, StatFit and MATWORKS [65] were used to obtain the adjustment and the parameter values of the density with reliability above 90 percent.

The corresponding density functions used are defined as follows [66]:

• Product: fresh chicken

  (a)

  Density function: gamma (G):

  $$f_X(x;\alpha ,\beta )=\frac{x^{\alpha -1}{e}^{-\beta x}}{\Gamma \left(\alpha \right){\beta }^{-\alpha }},0<x\le \infty ,\alpha =20,\beta =0.8$$

  (32)

  (b)

  Cumulative distribution function:

  $$F_X(u;\alpha ,\beta )={\int }_0^{x}f(u;\alpha ,\beta )du=\frac{\gamma (\alpha ,\beta x)}{\Gamma \left(\alpha \right)}=1-\sum _{i=0}^{\alpha -1}\frac{1}{i!}{\left(\frac{x}{\beta }\right)}^i{e}^{-\frac{x}{\beta }}$$

  where $\frac{\gamma (k,\frac{x}{\theta })}{\Gamma \left(\alpha \right)}$ is the lower incomplete gamma function.

  (c)

  Mean and variance:

  $${\mu }_X=\frac{\alpha }{\beta },{\sigma }^2=\frac{\alpha }{{\beta }^2}$$

  (d)

  Method to draw lotteries (using the Cheng’s algorithm) [67]:

  i.  

  Sample $U_1$ and $U_2$ from $\mathcal{U}$(0,1)

  ii. 

  $V\leftarrow \lambda ln\left[U_1/(1-U_2)\right]$

  iii.

  $X\leftarrow \alpha e^V$

  iv.

  if $b+d-X\ge ln\left[U_1^2{U}_2\right]$ deliver $X=\xi$

  v.  

  Go to step 1

  where $b=\alpha -ln\left(4\right)$, $d=\alpha +\frac{1}{\lambda }$ and $\lambda ={(2\alpha -1)}^{1/2}$

• Product: fresh fish

  (a)

  Density function: Weibull (W):

  $$f_X(x;\alpha ,\beta )=\frac{\alpha }{{\beta }^{\alpha }}{x}^{\alpha -1}{e}^{-{\left(\frac{x}{\beta }\right)}^{\alpha }},x\ge 0,\alpha =2.5,\beta =40$$

  (33)

  (b)

  Cumulative distribution function:

  $$F_X(u;\alpha ,\beta )=1-e^{-{\left(\frac{x}{\alpha }\right)}^{\beta }},x\ge 0$$

  (c)

  Mean and variance:

  $${\mu }_X=\beta \Gamma \left(1+\frac{1}{\alpha }\right),{\sigma }_X^2={\beta }^2\left[\Gamma \left(1+\frac{2}{\alpha }\right)-{\left(\Gamma \left(1+\frac{1}{\alpha }\right)\right)}^2\right]$$

  (d)

  Method to draw lotteries: inverse transformation:

  $$X=\xi =\beta {(-lnU)}^{1/\alpha },U\sim \mathcal{U}(0,1)$$

  (34)

• Product: fresh beef

  (a)

  Density function: Rayleigh (R):

  $$f_X(x;\sigma )=\frac{x}{{\sigma }^2}{e}^{-\frac{x^2}{2{\sigma }^2}},x\ge 0,{\sigma }^2=196$$

  (35)

  (b)

  Cumulative distribution function:

  $$F_X(x;\sigma )=1-e^{-\frac{x^2}{2{\sigma }^2}}x\ge 0$$

  (c)

  Mean and variance:

  $${\mu }_X=\sigma \sqrt{\frac{\pi }{2}}=1.253\sigma ,{\sigma }_X^2=\left(2-\frac{\pi }{2}\right){\sigma }^2=0.429{\sigma }^2$$

  (d)

  Method to draw lotteries: inverse transformation:

  $$X=\xi =\sigma \sqrt{-2lnU},U\sim \mathcal{U}(0,1)$$

  (36)

• Raw sausages, turkey and pork

  (a)

  Density function: normal (N):

  $$f_X(x;\mu ,{\sigma }^2)=\frac{1}{\sigma \sqrt{2\pi }},x\in \mathbb{R},\mu =19,\sigma =5$$

  (37)

  (b)

  Cumulative distribution function:

  $$F_X(x;\mu ,{\sigma }^2)=\frac{1}{2}\left[1+\mathrm{erf}\left(\frac{x-\mu }{\sigma \sqrt{2}}\right)\right],x\in \mathbb{R}$$

  (c)

  Mean and variance:

  $${\mu }_X=\mu ,{\sigma }_X^2=\sigma$$

  (d)

  Method to draw lotteries: simple composition:

  $$X=\xi =\mu +\sigma \left(\sum _{i=1}^{12}{U}_i-6\right),U_i\sim \mathcal{U}(0,1)$$

  (38)

• Fresh vegetables and eggs

  (a)

  Density function: Laplace (L):

  $$f_X(x;\mu ,b)=\frac{1}{2b}\exp \left(-\frac{|x-\mu |}{b}\right),\mu =40,b=4,t\ge 0$$

  (39)

  (b)

  Cumulative distribution function:

  $$F_X(x;\mu ,b)=\left\{\begin{array}{cc}\frac{1}{2}\exp \left(-\frac{\mu -x}{b}\right),\hfill & x<\mu \hfill \\ 1-\frac{1}{2}\exp \left(-\frac{x-\mu }{b}\right),\hfill & x\ge \mu \hfill \end{array}\right.$$

  $$=0.5[1+\mathrm{sgn}(x-\mu \left)\right(1-\exp (-|x-\mu |/b)\left)\right]$$

  (c)

  Mean and variance:

  $${\mu }_X=\mu ,{\sigma }_X=2b^2$$

  (d)

  Method to draw lotteries: inverse transform method:

  $$X=\xi =\mu -b\left[\mathrm{sgn}\left(U\right)ln(1-2|U\left|\right)\right],U\sim \mathcal{U}(0,1)$$

  (40)

• Itinerary on primary and secondary routes

  (a)

  Density function: exponential (E):

  $$f_X(x;\lambda )=\lambda e^{-\lambda t},t\ge 0$$

  (41)

  (b)

  Cumulative distribution function:

  $$F_X(t;\lambda )=1-e^{-\lambda t},t\ge 0$$

  (c)

  Mean and variance:

  $${\mu }_X=\frac{1}{\lambda },{\sigma }_X^2=\frac{1}{{\lambda }^2}$$

  (d)

  Method to draw lotteries: inverse transform method:

  $$X=\xi =-\frac{1}{\lambda }ln\left(U\right),U\sim \mathcal{U}(0,1)$$

  (42)

4. Numerical Example

Random search (RS) algorithms constitute a set of numerical optimization methods that do not require the use of gradient vectors as a sufficient condition to achieve convergence and approximate solutions to stochastic mathematical programming models, so they can be used in functions whereby they are not continuous or differentiable.

These classes of methods constitute the first attempts to use heuristics based on computer power to obtain approximate values for complex problems by using simple methods of progressing with a certain order or pattern in the parameter search space. The algorithm developed in our work works based on local random search, where each iteration depends on the candidate solution of the previous iteration. However, we have implemented a random search process that samples the entire search space (pure random search or global uniform random search).

The progression towards a solution is based on the creation of scenarios obtained by drawing all the random variables that intervene in it, evaluating at each step the value of the objective function and discarding the points that present feasibility but are not worthy to be incorporated because they do not contribute to the process of the maximization of the utility function. The convergence of the method and its efficiency are discussed in the appendix of this document.

In practice, there is a large number of probabilistic methods applied to the analysis of flexible manufacturing systems that try to reflect a more credible reality by incorporating uncertainty, see for example [68,69]. In particular, the design of distribution networks for PPs (green cold chain) presents a great challenge when uncertainty is included due to the large number of variables involved and their own nature. In our approach, we develop and apply a Generic Random Search Algorithm (Algorithm 1) to achieve an approximate solution to the proposed instance.

The Algorithm and Information Used

The algorithm works by using a sequence of iterations $\left\{x_k\right\}$ where in the iteration $k=0,1,\dots$, the new value of $x_k$ depends on the previous points and the parameters of the model. Therefore, in each k-th iteration, we have a feasible solution x and gradually progress toward a new solution. The Generic Random Search Algorithm used is shown below. Here, S is a feasible set for Equations (20)–(24), and for $\zeta =0$, we must have the real components of the following vectors:

$${\Theta }_0=(P,S,P+S,G,W,R,N,L,E,D,V),x_0=\left(t_0^{\prime }\right).$$

where all variables are defined below.

Algorithm 1 Generic Random Search Algorithm

Input: Initialize the parameters ${\Theta }_{\zeta }$ of the algorithm with the initial set of points $x_0=$  $\{x_{\zeta }^1,x_{\zeta }^2,\dots ,x_{\zeta }^k\}$ for each $x_{\zeta }^i\subset S$, and iteration index $\zeta =0$.

Output: $Z_0^{*}\left(x_0^{*}\right)=max\{Z_0\left(x_{\zeta }^i\right):x_{\zeta }^i\in x_0\}$  $\zeta \leftarrow 1$  $x_{\zeta }\leftarrow x_0^{*}$  while $ϵ\le 0.01$ do    Generate a $x_{\zeta }=\{x_{\zeta }^1,x_{\zeta }^2,\dots ,x_{\zeta }^k\}$ for each $x_{\zeta }^i\subset S$    $Z_{\zeta }^{*}\left(x_{\zeta }^{*}\right)\leftarrow max\{Z_{\zeta }\left(x_{\zeta }^i\right),Z\left(x_{\zeta }\right)\mid x_{\zeta }^i\in x_{\zeta }\}$   if $Z_{\zeta }^{*}\left(x_{\zeta }^{*}\right)<Z_{\zeta }\left(x_{\zeta }^i\right)$ then $x_{\zeta +1}\leftarrow x_{\zeta }^{*}$,   else if $Z_{\zeta }^{*}\left(x_{\zeta }^{*}\right)\ge Z_{\zeta }\left(x_{\zeta }^i\right)$ then $x_{\zeta +1}\leftarrow x_{\zeta }^i$   end if   $ϵ\leftarrow E\mid Z\left(x_{\zeta +1}\right)-Z\left(x_{\zeta }\right)\mid$   $\zeta \leftarrow \zeta +1$   Update the ${\Theta }_{\zeta }$ parameter values  end while

where S is a feasible set for Equations (20)–(24), and for $\zeta =0$ we must have the real components of the following vectors:

$${\Theta }_0=(P,S,P+S,G,W,R,N,L,E,D,V),x_0=\left(t_0^{\prime }\right)$$

with components given by

• P: primary path travel time (in hours).
• S: secondary path travel time (in hours).
• G: sample value draw for the gamma distribution.
• W: sample value draw for the Weibull distribution.
• R: sample value draw for the Rayleigh distribution.
• N: sample value draw for the normal distribution.
• L: sample value draw for the Laplace distribution.
• E: sample value draw for the exponential distributions.
• D: sample value draw for the consolidated demand.
• V: number of vehicles required with a capacity of 30 tons to transport demand D.
• $t^{\prime }$: time window required for shipment.
• Z: objective function value.

The consecutive values of the sequence are selected in such a way that

$$Z\left(x_1\right)\le Z\left(x_2\right)\le \dots \le Z\left(x_{\zeta }\right),$$

and at the optimal value, it is verified that

$$Z\left(x^{*}\right)\ge Z\left(x_{\zeta }\right),\forall \zeta$$

The update procedure is as follows:

$$x_{\zeta +1}=\left\{\begin{array}{cc}{x}_{\zeta +1}^{\prime }\hfill & \mathrm{if}Z\left(x_{\zeta +1}^{\prime }\right)<Z\left(x_{\zeta }\right)\hfill \\ x_{\zeta }\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(43)

Typically, for the use of random search algorithms for the stochastic optimization of models of the type shown in Equation (1), the following stopping rules are common:

• For a given number of iterations $\zeta$,

  $$minE(\parallel x_{\zeta }-x^{*}\parallel )=E(\mid x_{\zeta }^k-x^{*}\mid )$$

  (44)
• When a predetermined $ϵ$-value is known,

  $$maxP(\parallel x_{\zeta }-x^{*}\parallel )\le ϵ=P(\parallel x_{\zeta }^k-x^{*}\parallel )\le ϵ$$

  (45)
• When a predetermined $ϵ$-value is known,

  $$maxP(\mid Z\left(x_{\zeta }\right)-Z\left(x^{*}\right)\mid )\le ϵ=P(\mid Z\left(x_{\zeta }^k\right)-Z\left(x^{*}\right)\mid )\le ϵ$$

  (46)
• Depending on the restrictive set

  $$min\zeta ={\zeta }^s,s=1,2,\dots S$$

  (47)

  subject to

  $$E\left(x_{\zeta }^s\right)\in [R_x=\{x:\parallel x-x^{*}\parallel \le {ϵ}_x\left\}\right]$$

  (48)

  or

  $$E\left(x_{\zeta }^s\right)\in [R_Z=\{x:\mid Z\left(x\right)-Z\left(x^{*}\right)\mid \le {ϵ}_Z\}]$$

  (49)

The stopping criterion that we use in our algorithm is given by the set of Equations (47) and (49). The aforementioned results were achieved before reaching 100 iterations in the process with a proposed value of ${ϵ}_Z=0.01$. In Appendix A, we discuss the efficiency of the proposed algorithm.

For the developed instance, we have the following information.

Table 2. Cumulative distribution functions $F_t\left(t^{\prime }\right)$ and the $P(W_t\le t^{\prime })$ function.

${\mathit{t}}^{\prime }$	G	W	R	N	L	$1-\left[{\Pi }_{\mathit{i}=1}^5[1-{\mathit{F}}_{{\mathit{T}}_{\mathit{i}}}\left({\mathit{t}}^{\prime }\right)]\right]$
2	0.0000	0.0000	0.0000	0.0003	0.0000	0.0003
4	0.0000	0.0031	0.0198	0.0013	0.0000	0.0242
6	0.0000	0.0086	0.0440	0.0046	0.0001	0.0568
8	0.0000	0.0177	0.0768	0.0139	0.0001	0.1060
10	0.0000	0.0307	0.1174	0.0359	0.0000	0.1756
12	0.0000	0.0481	0.1647	0.0807	0.0000	0.2694
14	0.0000	0.0699	0.2172	0.1586	0.0000	0.3879
16	0.0009	0.0962	0.2738	0.2742	0.0012	0.5247
18	0.0034	0.1270	0.3330	0.4207	0.0020	0.6645
20	0.0102	0.1620	0.3934	0.5792	0.0033	0.7890
22	0.0251	0.2009	0.4539	0.7257	0.0055	0.8839
24	0.0524	0.2432	0.5132	0.8413	0.0091	0.9451
26	0.0961	0.2886	0.5704	0.9192	0.0151	0.9780
28	0.1580	0.3363	0.6246	0.9640	0.0248	0.9926
30	0.2371	0.3856	0.6753	0.9860	0.0410	0.9979
32	0.3295	0.4358	0.7219	0.9953	0.0776	0.9995
34	0.4292	0.43863	0.7642	0.9986	0.1115	0.9999
36	0.5297	0.5363	0.8020	0.9996	0.1839	0.9999
38	0.6248	0.5850	0.8355	0.9999	0.3032	0.9999
40	0.7100	0.6321	0.8646	0.9999	0.5000	0.9999
42	0.7826	0.6768	0.8897	0.9999	0.6967	1.0000
44	0.8417	0.7189	0.9110	0.9999	0.8160	1.0000
46	0.8880	0.7578	0.9289	0.9999	0.8884	1.0000
48	0.9227	0.7935	0.9438	0.9999	0.9323	1.0000
50	0.9480	0.8256	0.9560	0.9999	0.9589	1.0000
52	0.9659	0.8544	0.9659	1.0000	0.9751	1.0000
54	0.9781	0.8797	0.9739	1.0000	0.9849	1.0000
56	0.9863	0.9016	0.9802	1.0000	0.9908	1.0000
58	0.9915	0.9205	0.9851	1.0000	0.9944	1.0000
60	0.9949	0.9364	0.9889	1.0000	0.9966	1.0000
62	0.9970	0.9498	0.9918	1.0000	0.9980	1.0000
64	0.9982	0.9608	0.9940	1.0000	0.9988	1.0000
66	0.9990	0.9697	0.9957	1.0000	0.9992	1.0000
68	1.0000	0.9769	0.9969	1.0000	0.9995	1.0000
70	1.0000	0.9826	0.9978	1.0000	0.9997	1.0000
72	1.0000	0.9871	1.0000	1.0000	1.0000	1.0000
74	1.0000	0.9905	1.0000	1.0000	1.0000	1.0000
76	1.0000	0.9931	1.0000	1.0000	1.0000	1.0000
78	1.0000	0.9951	1.0000	1.0000	1.0000	1.0000
80	1.0000	0.9965	1.0000	1.0000	1.0000	1.0000
82	1.0000	1.0000	1.0000	1.0000	1.0000	1.0000

shows the numerical results of the distributions involved in the process for $\theta =5$ and $2\le t^{\prime }\le 82$ (in hours), with the parameter values indicated in Equations (32), (33), (35), (37) and (39).

The numerical values used for $\lambda$ in the primary and secondary tours are described in

Table 3. Numerical values used for $\lambda$ in the primary tours and loading + unloading times.

	${\mathcal{T}}_0$	${\mathcal{T}}_1$	${\mathcal{T}}_2$	${\mathcal{T}}_{12}$	${\mathcal{T}}_{23}$	${\mathcal{T}}_{34}$	${\mathcal{T}}_{45}$
$\lambda$	3	2.5	3	3.2	3	3.5	3
$\mathfrak{d}$ (in hours)				1.5	1	1.2	1.3

and

Table 4. Numerical values used for $\lambda$ in the secondary tours and loading + unloading times.

							$\mathfrak{d}$
							(in hours)
	${\tau }_{12}^1$	${\tau }_{23}^1$	${\tau }_{34}^1$	${\tau }_{45}^1$	${\tau }_{56}^1$	${\tau }_{67}^1$
$\lambda$	1.6	2.1	3.2	2.1	2.6	3.5	1.2
	${\tau }_{12}^2$	${\tau }_{23}^2$	${\tau }_{34}^2$	${\tau }_{45}^2$	${\tau }_{56}^2$	${\tau }_{67}^2$
$\lambda$	1.5.	2.1	2.3	1.3	1.6	2	1.2
	${\tau }_{12}^3$	${\tau }_{23}^3$	${\tau }_{34}^3$
$\lambda$	2.3	2.5	2.9				0.8
	${\tau }_{12}^4$	${\tau }_{23}^4$	${\tau }_{34}^4$	${\tau }_{45}^4$	${\tau }_{56}^4$	${\tau }_{67}^4$
$\lambda$	1.2	2.5	2.5	1.9	1.6	1.7	0.5

. The corresponding demand values in the shoveling nodes are generated by using the procedure described in Equation (38).

Table 5. Numerical values used for the draw of demands.

Shoveling Node	Mean ($\mathit{\mu }$)	Variance $\left({\mathit{\sigma }}^2\right)$
$d_1$	30	10
$d_2$	60	14
$d_3$	20	9
$d_4$	100	18
Total	210	51

shows the parametric values of the distribution used to draw the demand values.

The associated information of the process is as follows. Table 1 shows the numerical results of the distributions involved in the process for $\theta =5$ and $2\le t^{\prime }\le 82$ (in hours), with the parameter values indicated in Equations (32), (33), (35), (37) and (39).

The numerical values used for $\lambda$ in the primary and secondary tours are described in Table 2 and Table 3. The corresponding demand values in the shoveling nodes are generated by using the procedure described in Equation (38). Table 4 shows the parametric values of the distribution used to draw the demand values.

By virtue of Equation (17), $D\sim \mathcal{N}\left({\sum }_i{\mu }_i,\sqrt{{\sum }_i{\sigma }^2}\right)$,

$$t^{\prime }=min\{T_s:max\left\{{\omega }_k\right\}\le min\left\{T_s\right\}\}$$

(50)

As an example, for $t^{\prime }=21.68$, we have $max\left\{{\omega }_k\right\}=max\{17.73,18.06,13.36,18.82\}=18.82$ and $min\left\{T_s\right\}=min\{21.68,46.68,42.54,26.42,35.09\}$ = 21.68. The value of Z is obtained by interpolating the value obtained from $t^{\prime }$ in Table 1. Equation (50) guarantees that the journey through the network will be conducted in a time less than or equal to the maximum travel time throughout the network. Finally, columns twelve and thirteen show the required value of the consolidated demand (with $\alpha =0.05)$ and the total number of vehicles necessary for its transportation. Finally,

Table 6. Results of 10 computational runs associated with the instance.

k	P	S	Total	G	W	R	N	L	${\mathit{t}}^{\prime }$	Z	D	V
1	3.80	11.68	15.48	19.84	42.71	73.84	27.98	34.46	19.84	0.7826	189	6
2	6.46	12.23	18.70
3	10.08	6.38	16.46
4	12.60	6.91	19.53
1	4.30	12.37	16.67	20.11	43.41	53.55	26.33	29.38	20.11	0.7933	216	7
2	8.72	8.23	16.95
3	10.14	4.74	14.88
4	13.89	5.50	19.38
1	2.59	5.53	8.12	20.23	46.62	56.67	26.52	31.57	20.23	0.7980	205	7
2	3.74	8.81	12.55
3	7.04	7.18	14.23
4	9.43	9.88	19.31
1	5.95	14.00	19.95	20.39	50.89	71.33	27.80	36.07	20.39	0.8034	184	6
2	7.27	11.72	18.99
3	8.91	4.86	13.78
4	10.70	7.07	187.77
1	5.22	13.88	19.11	20.61	48.10	84.55	24.83	31.39	20.61	0.8131	248	8
2	6.61		12.18
3	9.90	5.53	15.43
4	11.34	7.26	18.61
1	3.79	13.29	17.08	20.65	53.62	60.88	28.29	30.22	20.65	0.8146	214	7
2	5.61	10.09	15.70
3	7.55	4.62	12.17
4	10.33	8.10	18.44
1	3.25	16.06	19.31	20.69	40.77	59.84	23.85	31.61	20.69	0.8162	215	7
2	4.37	8.02	12.40
3	6.49	4.81	11.31
4	9.18	9.8	18.98
1	1.96	15.89	17.85	20.92	47.81	61.33	25.54	32.18	20.92	0.8252	194	7
2	9.26	7.47	16.73
3	11.47	5.52	17.00
4	13.46	6.41	19.88
1	4.97	10.90	15.88	20.98	44.85	56.79	28.43	29.36	20.98	0.8276	214	7
2	6.74	14.12	20.86
3	11.00	5.32	16.32
4	13.66	6.04	19.71
1	2.71	15.01	17.73	21.68	46.68	42.54	26.42	35.09	21.68	0.8347	194	7
2	6.80	11.25	18.06
3	8.14	5.21	13.36
4	15.90	2.91	18.82

shows the 10 convergence results of the algorithm.

Figure 3 shows the convergence toward a local solution of the algorithm. Thus, for this instance, the optimal transportation plan is given by achieving a route with a time window less than or equal to 21.68 h, hoping to achieve it with a probability of 0.8347.

5. Conclusions

In this research, we addressed the real problem of designing strategies for the optimal management of distribution routes for perishable products that involve a time window. Our objective was to maximize the probability of meeting the time limit assigned to the routing problem with a fixed capacity, thus affecting the issue of minimizing the pollution generated by refrigerated transportation.

Our assumptions consider the existence of primary and secondary routes as well as spray nodes along the trajectories inherent to the distribution network. To give more reality to the model, we considered that the travel times in the network are random, as well as the demand and shelf life of the transported products. The resulting model is highly dynamic and stochastic, which is why we designed a solution strategy based on Monte Carlo optimization by a random search technique. The results obtained suggest that the design of networks like the one shown here or perhaps more complex ones are feasible to address by using this approach by adapting search algorithms (heuristics) that incorporate simple convergence rules.

The random search method used to obtain a local solution to the proposed instance can be considered as blind search [70] since the algorithm repeatedly samples the feasible region S by means of one or several theoretical or empirical probability distributions. Following this procedure, each candidate point is updated by generating a set of sample values of the distributions involved in the process in S and is updated only if the candidate point is improving. The convergence in probability (within an $ϵ$ neighborhood) toward a global optimum of a random search process was proved before and is considered as the global phase in clustering and multiple-start algorithms [71,72].

The limitations of our work are the following: (a) historical information on the shelf life of the products analyzed must be available in order to allow for statistical analysis and Monte Carlo simulation; (b) the previous issue also applies to the evaluation of travel times and vehicle loading and unloading operations, breaks, meals and more; (c) knowledge of historical demand information is essential to evaluate the quantities to be shipped; and (d) ignorance of information in this regard invalidates the assumptions of the model.

Future research should involve the computational complexity associated with more complicated topologies and other cases of interest such as transportation with cargo collection, volumetric of transportation units, reverse collection and limitations on the availability of transportation and operators.