BI-OBJECTIVE APPROACHES TO DEAL WITH ACCIDENT RISK AND LOGISTIC COSTS IN VEHICLE ROUTING PROBLEMS

ABSTRACT

Vehicle Routing Problems (VRP) have been widely researched throughout history as a way of optimising routes by minimising distances and planning deliveries efficiently, but the issue of risk in VRP has received less attention over time. This is essential to increasing transport safety to avoid interruptions in supply chains and improve delivery reliability. Therefore, this study aims to support decision makers to plan routes for road freight transportation considering not only the logistics cost, but also travel safety due to road hazards. An analytical approach based on statistics was developed in which data of vehicle accidents and road features were used to estimate the risk cost by using the Monte Carlo simulation. A bi-objective approach based on PROMETHEE II and ε-constrained methods was used in the Capacitated Vehicle Routing Problem (CVRP) to analyse the conflict between the logistic cost and accident risk. Key contributions of this study are an analytical approach based on statistics, Monte Carlo simulation and multicriteria methods in a CVRP model to explore the trade-off between logistic costs and accident risk expressed by a risk cost. The outcomes of this study show to be useful in practice to analyse transportation decisions in the VRP model involving route planning considering accident risk.

Keywords:
CVRP; VRP with risk; Monte Carlo simulation; Promethee II; ε-constrained

1 INTRODUCTION

In an increasingly globalised and connected world in which the population and urbanisation increase every year, the pressure for more efficient, sustainable and safe road freight has risen due to the greater demand for different types of products. Road freight plays a fundamental role in global logistics due to the flexibility of the infrastructure that this segment of transport offers, allowing a door-to-door service, which is not possible in maritime, air and rail transportation modes (EngstrÖm, 2016).

On the other hand, the reliability in road freight depends on some elements that can be the cause of many accidents, interrupting the supply chain and generating significant losses. According to data from the World Health Organization (WHO), road accidents cost on average around 3% of the Gross Domestic Product (GDP) of a country and lead to the death of 1.35 million people, and is the eighth leading cause of death all around the world.

Thus, reducing the accident rate of road freight is an important factor from a strategic and sustainable point of view, which aims to ensure logistical and economical development. When an accident occurs, financial losses can lead to large financial proportions, in addition to causing interruptions in the supply chain. From this point, the need to reduce losses in road transportation can be highlighted (EngstrÖm, 2016).

When a route is chosen for a vehicle that will leave a depot to deliver goods to other locations, not only the logistic costs and the distance should be considered, but also some route hazard measures should be added. Most of the VRP studies consider only the logistic cost in the models to solve the routing problem. Recently risk measurement was addressed in cash-in-transit problems through the VRP, which consisted of a model that aimed to minimise the distances travelled with the restriction that the risk value of robberies of heavy trucks during the transportation of money was limited by a risk threshold (Talarico et al., 2015). More recently the route safety was also discussed in studies on hazardous materials transportation such as fuel, flammable materials, gases and others, aiming to reduce social and environmental impacts and increase transport safety by minimising the risk factor (Holeczek, 2021).

Despite the fact that the VRP has been widely studied throughout its sixty-year history, few studies have addressed different types of risk in the problem. As mentioned, some research has emerged in the areas of cash-in-transit and hazardous materials, but they did not consider accident risk due to road conditions in the models. This is likely because in most cases real data is difficult to find and a methodology to estimate the risks should be developed.

Therefore, verifying the applications of the subject risk in VRP and considering the lack of studies that address the accident risk due to characteristics of the routes (infrastructure, traffic etc.), this study intends to contribute to filling this gap by developing an analytical approach based on simple statistical calculations coupled with two multi-criteria methods in a CVRP model to analyse the trade-off between the logistics cost and accident risk measured by cost. The Preference Ranking Organization Method for Enrichment Evaluations (PROMETHEE) II and the ε-constrained methods were selected to solve this conflicting bi-objective problem. Finally, the statistical approach to estimate risk cost was validated with experts in cargo insurance and the obtained solutions were interpreted knowing the characteristics of the roads well in the tested instance, indicating that the outcomes were compatible with those expected.

This paper is organized in four sections, besides this introduction. Section 2 reviews and discusses key studies which address risks in VRP. In Section 3, the proposed analytical approach is presented in detail. Finally, the results (Section 4) are presented through the experimental studies and the conclusion is drawn in Section 5.

2 LITERATURE RELATED TO RISK IN VRP

Vehicle Routing Problems (VRP) have been extensively studied throughout their history to support real-life applications. Risk and travel safety in VRP have received greater attention in applications for road freight of hazardous materials, whose risk is associated with socio-environmental damage, and cash-in-transit, which is related to cargo theft (Talarico et al., 2017).

Erkut & Ingolfsson (2005) cited eight risk models that were developed for the optimisation of hazardous materials transportation, whose three main ones are presented in Table 1. Considering a path r as a set of links {i ₁ , i ₂ , ..., i _n }, the first risk model IP _r is calculated through the probability of the undesirable occurrence event of each link p _i , while the second PE _r considers only the number of people exposed to risk D _i . In the traditional model, risk TR _r is a product between the probability of the undesirable event p _i and its measure of consequence C _i .

Table 1
Risk assessment models. Adapted from Erkut & Ingolfsson (2005).

Model	Equation
Incident Probability	$I{P}_r=\sum _{i\in r}{p}_i$	pi = accident probability
Population Exposure	$P{E}_r=\sum _{i\in r}{D}_i$	Di = population exposure
Traditional Risk	$T{R}_r=\sum _{i\in r}{p}_i.C_i$	pi = accident probability
		Ci = measure of the consequence

To investigate the behaviour of these three risk models in VRP, Holeczek (2021) used monoobjective functions that minimise distance, accident risk and population exposure.

The results showed that the Traditional Risk TR _r generated the lowest total risk value when compared to the other models, however when considering the total distance obtained for the TR _r , the greatest deviation in relation to the minimum total distance can be observed.

The Incident Probability IP _r offers the best trade-off with an economical goal and it is most appropriate for problems where the consequences are uncertain. Regarding the Population Exposed PE _r , the data are more easily acquired and the results are evaluated more intuitively by the decision maker, but it can only be applied to problems that consider urban areas, because for an environment, such as rural areas, other factors must be considered.

The elements of each model such as accident probability, consequences, number of people exposed and others, are defined according to the case being studied. Androutsopoulos & Zografos (2012) and Carrese et al. (2022), for example, studied risks in hazardous materials transportation and they considered that the undesirable event would be the road accident whose consequence is related to the number of people exposed to risk. On the other hand, Talarico et al. (2015) applied the risk in the cash-in-transit routing problem that they considered the probability of a robbery is proportional to the distance of the route and, as a consequence, the amount of cash transported in the truck.

In addition to those elements, the models also vary according to the risk approach used in VRP. Thus, Table 2 shows some key features in the risk calculations used by the authors.

Du et al. (2017), Pradhananga et al. (2014) and Wang et al. (2018) applied the concept of the traditional model, as proposed by Erkut & Ingolfsson (2005), in which the probability of an accident was used as an undesirable event and the exposed population as a consequence. Androutsopoulos & Zografos (2012) and Holeczek (2021) applied the traditional definition, however it was considered as load dependent, that is, the amount of load factor is added to the model and varies as deliveries are made. Carrese et al. (2022) was also addressed the traditional risk model, but two other factors that interfere in the driver’s attention are added to the objective function: the Altimetric Index and the Planimetric Index. The first considers the ground elevations along the route while the second is introduced to take into account geometrical constraints related to the road radius.

Bula et al. (2016, 2019) also uses the traditional risk which comprises the probability of an undesirable event p _i and population exposure as the measure of consequence C _i . They consider several aspects that play an important role in determining the likelihood of an undesirable occurrence. Therefore, the p _i is combined as a result of: accident probability related to the type of truck, probability of hazardous materials released in case of an accident, parameters that represent the characteristics of the materials and the amount of cargo carried, and also the length of arc in the route.

Chai et al. (2023) argue that many studies consider only issues such as accident probabilities, the population exposed and consequences for the environment, but it is also important to consider the driver’s behavior as a main cause of accidents. Thereby, Chai et al. (2023) introduce the driver’s behavior as a factor that also influences the risk of accidents in the transport of hazardous materials. Some aspects that impact the driver risk were considered, such as: age, driving experience, educational background, gender, driving speed and driving habits.

As a variation of hazardous materials models, the cash-in-transit ones arise. Talarico et al. (2015, 2017) proposed the traditional method to calculate the risk of robbery. As already mentioned, they considered the consequence as equivalent to the amount of value being transported. Ghannadpour & Zandiyeh (2020) also assumed the consequence as the same manner and the distance travelled is proportional to the risk of theft, but added to the model a factor relative to the frequency of using the same route and the ambushing probabilities to the vehicle and its success. The authors also state that the probability of a robber attack was estimated using game theory and a model minimizes the risk of cash-in-transit developed using multi-criteria decision-making approaches.

Talarico et al. (2015) and Bula et al. (2016) used a mono-objective function in which the first authors minimised the distances and the second the accident risk. Androutsopoulos & Zografos (2012), Bula et al. (2019), Pradhananga et al. (2014), Ghannadpour & Zandiyeh (2020) and Wang et al. (2018) suggested bi-objective functions that analyse both logistics costs or distances and risks. Multi-objectives are presented by Carrese et al. (2022) and Zheng (2010) whereby the last one aimed to minimise the distance, weighted the risk by the traditional method (probability of accident and a consequence) adding the number of people exposed.

Branko Milovanović (2012) developed a methodology to calculate road accident risks when transporting hazardous materials, which considers some elements that influence the accident probability, as well as elements that influence their consequences. These elements were measured from indirect interviews in which experts obtained numerical risk results for each route through this analysis. On the other hand, it did not use mathematical models of VRP to optimise routes and it did not consider statistical analysis.

Generally, the risk factor is analysed as an objective to be minimised, however Talarico et al. (2015) considered it as a constraint in which the risk value is defined by a risk threshold and classified as the Risk constrained Cash-in-Transit Vehicle Routing Problem (RCTVRP). Wang et al. (2018) restricted the condition that no vehicles of the same fleet should travel in echelon because when there are two or more vehicles using the same route at the same time, the consequences are considered to be greater whether or not an accident occurs between them.

For risk measurement, some studies addressed how the data were explored, and according to Du et al. (2017), real historical data from accidents should be integrated and resort to big data in terms of formulating models to transport hazardous materials. Moreover, Talarico et al. (2015) mention that there is lack of data available and Androutsopoulos & Zografos (2012) indicate there is no data exploration related to risk measurement due to its complexity and also states that future studies should deal with this issue.

Pradhananga et al. (2014) estimated accident rates using data collected from the Institute for Traffic Accident Research and Data Analysis (ITARDA) and the Ministry of Land, Infrastructure, Transport and Tourism (MLIT), both from Japan. For a future study, Pradhananga et al. (2014) proposed extensions of the model considering such characteristics for the hazardous materials routing problem using real-time traffic information and the effects of infrastructural characteristics of the road network.

Carrese et al. (2022) calculated the road accident probability using data obtained by the mobility agency in Rome, quantified population density through a census data and measured the infrastructure through the Google Application Programming Interface (API).

Although some studies still try to deal with using real data in their problems, the probabilistic view in data analysis is not addressed in more depth. This study aims to use the probabilistic approach of the accident occurrences in the model. Table 2 summarises the problem features such as the VRP type, objective function and real world-data with a VRP taxonomic discussed in Braekers et al. (2016). It can be observed that the studies are concentrated in only two areas: hazardous materials and cash-in-transit. Regarding other segments of transport, there are few studies that take into account risks in VRP despite its importance as a result of damage when an accident occurs.

Table 2
Characteristics of risks in VRP studies.

Authors	Case Studied	VRP type	Objective Function	Risk Approach	Real data
Androutsopoulos & Zografos (2012)	Hazardous Materials	Vehicle Routing Problem with Time Windows	Travel time dependent and Implied hazard/risk related	accident probability * population exposure * load amount	No
Bula et al. (2016)	Hazardous Materials	Homogeneous Vehicle Routing Problem	Implied hazard/risk related	accident probability related to truck type * release probability * load characteristics * route length * load amount population exposure	No
Bula et al. (2019)	Hazardous Materials	Homogeneous Vehicle Routing Problem	Distance dependent and Implied hazard/risk related	accident probability related to truck type * release probability * route length * load characteristics * load amount * population exposure	No
Carrese et al. (2022)	Hazardous Materials	Vehicle Routing Problem with Time Windows	Travel time dependent and Implied hazard/risk related	accident probability * population exposure + altimetric index + planimetric index	Yes
Du et al. (2017)	Hazardous Materials	Multi-depot Vehicle Routing Problem	Implied hazard/risk related	accident probability * population exposure	Yes
Ghannadpour & Zandiyeh (2020)	Cash-in-Transit	Vehicle Routing Problem with Time Windows	Distance dependent and Implied hazard/risk related	robbery attack probability * theft success probability * route length * load amount * frequency of repeated use of a route	No
Holeczek (2021)	Hazardous Materials	Capacitated Vehicle Routing Problem	Distance dependent and Implied hazard/risk related	accident probability * population exposure * load amount	Yes
Pradhananga et al. (2014)	Hazardous Materials	Vehicle Routing Problem with Time Windows	Travel time dependent and Implied hazard/risk related	accident probability * population exposure	Yes
Talarico et al. (2015)	Cash-in-Transit	Risk constrained Cash-in-Transit Vehicle Routing Problem	Distance dependent	arc length * load amount	No
Talarico et al. (2017)	Cash-in-Transit	Risk constrained Cash-in-Transit Vehicle Routing Problem	Distance dependent	arc length * load amount	No
Wang et al. (2018)	Hazardous Materials	Vehicle Routing Problem with Time Windows	Distance dependent and Implied hazard/risk related	accident probability * population exposure	No
Zheng (2010)	Hazardous Materials	Capacitated Vehicle Routing Problem	Distance dependent, Implied hazard/risk related and Others	accident probability * consequence + population exposure	No
Chai et al. (2023)	Hazardous Materials	Vehicle Routing Problem with Soft Time Window for Hazardous Materials	Distance dependent and Implied hazard/risk related	accident probability * population exposure * driving risk * load amount * arc length	Yes
This study	Indistinct load type	Capacitated Vehicle Routing Problem	Distance dependent and Implied hazard/risk related	accident probability, road infrastructure and traffic, load value	Yes

3 ANALYTICAL APPROACH

In this section, a detailed description with an illustration of the workflow to apply the analytical approach is provided (Figure 1). Firstly, the parameters, variables and constraints of the model representing the Capacitated Vehicle Routing Problem (CVRP) were defined and the logistic costs (c _{i j} ) were calculated by the online tool Qualp. (2021), which considers fuel expenses, based on the vehicle’s consumption, and tolls along the road if applicable.

The accident probability (Paccident _{i j} ) was generated by a statistical calculation that considers historical records from Brazilian government agencies whose information is publicly available on the web and from a load insurance company. Among the government agencies, there are: the National Transport Confederation (CNT), the Department of Roads and Highways of the State of Sao Paulo (DER-SP) and the National Department of Transport Infrastructure (DNIT). Some characteristics that interfere in the accident probabilities were extracted by Branko Milovanović (2012) and all collected and processed data used in this study is available at Bilato (2022).

The data were processed using the Knime Analytic Platform tool, which generated the accident probabilities for each arc. Then, from the probability results, the Monte Carlo simulation was implemented to obtain the risk costs (r _{i j} ).

Two different methods were performed to deal with the bi-objective problem, i.e. the Preference Ranking Optimization Method for Enrichment Evaluation (PROMETHEE II) and the εconstrained method. The CVRP and these two methods were implemented using Python, and the results were obtained by Gurobi Optimization.

The results from PROMETHEE II and ε-constrained methods were compared and analysed. The method used to calculate the r _{i j} will be explained and how the PROMETHEE II and εconstrained method were used in the CVRP.

Calculating the accident risk cost

As already mentioned, the risk cost r _{i j} was generated for each arc connecting location i to location j. A location can be a city or centroid, for example. The arc ij may have more than one distinguished roads.

These calculations were performed by using the Knime Analytics Platform, according to the workflows represented in Figures 2, 3 and 4 and Equations (1), (2), (3), (4), (5), (6), (7), (8) and (9). This approach to estimate r _{i j} was necessary particularly because there are little data available referring to road accidents.

A general probability of accidents (Pgeneral) occurring on any road in Brazil according to the workflow shown in Figure 2 was estimated. Data collected from DER-SP (2021) indicate the flow of all vehicles V _sp and heavy vehicles HV _sp in the State of Sao Paulo (sp). They were used to calculate the share of heavy vehicles P _sp that circulates in the roads of Sao Paulo State by Equation (1).

Figure 1
Workflow of the analytical approach in this study.

Figure 2
Workflow to calculate Pgeneral.

For the sake of simplicity and lack of specific data, it is assumed that P _sp is the same for all Brazilian Federal roads, therefore in Equation (2) the flow of heavy vehicles (HV ) is calculated using the flow of all vehicles in Federal roads (V _br ) - data extracted from DNIT (2021). Finally, Pgeneral was generated by Equation (3) from the number of accidents in Brazilians’ Federal roads (N _accidents ), collected from CNT (2019), and HV from Equation (2).

$P_{sp}=H{V}_{sp}/V_{sp}$ (1)

$HV=P_{sp}.V_{br}$ (2)

$Pgeneral=\frac{N_{accidents}}{HV}.100$ (3)

The elements considered for the calculation of r _{i j} are: the types of roads and the traffic of heavy vehicles. These elements were based on Branko Milovanović (2012) and the first is considered because in Brazil there are several types of roads that demonstrate different safety levels. The Accident Panel Report prepared by CNT (2019) breaks it into five categories as presented by Table 3. The number of deaths per hundred of accidents extracted from this report is also shown, and is used in the calculations for accident risk.

Table 3
Death rate per type of road CNT (2019).

Road type	Death rate per 100 accidents
Two-lane two-way road with central safety lane	12.3
Two-lane two-way road with central barrier	8.5
Two-lane two-way road with central line	18.0
Single-lane one-way road	11.9
Single-lane two-way road	22.3

The flow of heavy vehicles is obtained by speed radars installed along the roads (DER-SP (2021)) and it was considered as directly proportional to the accident probabilities as it was observed that two roads had the same infrastructure, distinguished only by the flow: the one with highest flow had more accidents, thus having greater probability of this occurrence.

Parameters it _h and iv _h represent the types of road (i.e. number of lanes, central barrier etc.) and the flow of heavy vehicles, respectively. These are calculated according to the workflow shown in Figure 3, where road h belongs to the problem roads of set H (h ∈ H) (CNT (2019) and DER-SP (2021)).

Initially, $\overline{x}$ and $\overline{y}$ must be calculated by Equations (4) and (5), which represent the average flow of vehicles x _h running in road h and the y _h the death rate of road type h. In the study, as there are 49 roads and five types of them, the parameter are Nh = 49 and Nt = 5.

Then, the parameters iv _h and it _h were calculated by Equations (6) and (7), which basically consists of a condition, i.e. if iv _h or it _h is greater than 1.0, the accident probability on road h will be greater than the general probability (P _general ), otherwise the values are less than 1.0.

It is important to explain that two or more arcs comprising the same roads may present different risk costs according to the road distance in those arcs. For example, in a given arc, the vehicle would be on the more dangerous road for longer than the other arc in which a vehicle would travel for a short time on this same dangerous road.

$\overline{x}=\frac{{\sum }_{h\in H}{x}_h}{Nh}$ (4)

$\overline{y}=\frac{{\sum }_{t\in T}{y}_t}{Nt}$ (5)

$i{v}_h=1+\frac{x_h-\overline{x}}{\overline{x}},\forall h\in H$ (6)

$i{t}_h=1+\frac{y_h-\overline{y}}{\overline{y}},\forall h\in H$ (7)

$e_{ij}=\frac{{\sum }_{h\in H}i{v}_h.i{t}_h.l_h}{l_{ij}},\forall (i,j)\in E$ (8)

$Pacciden{t}_{ij}=Pgeneral.e_{ij},\forall (i,j)\in E$ (9)

Figure 3
Workflow to calculate parameters ivh e ith.

In some arcs, the vehicle may pass through more than one road (h), i.e. the arc Limeira to Cosmópolis has two roads: h ₁ named SP330 and h ₂ named SP133 whose flows and features are distinguished. Thus, a parameter named e _{i j} is needed, which is weighted by iv _h , it _h and l _h , and represented by Equation (8) as shown in the workflow of Figure 4. Where l _h is the length that the truck traverses each road h of the arc i j and l _{i j} represents the total length of the arc. Finally, Equation (9) describes the accident probability Paccident _{i j} of the arc i j.

Figure 4
Workflow to calculate ei j and Paccidenti j .

Table 4
Value involved in accident cargo transportation.

Range of values	occurrence
$ 0.01 to $ 200,000.00	37.91%
$ 200,000.00 to $ 300,000.00	24.17%
$ 300,000.00 to $ 500,000.00	19.91%
$ 500,000.00 to $ 1,000,000.00	16.11%
$ 1,000,000.00 or more	1.90%

After finding the Paccident _{i j} , it is possible to estimate r _{i j} using the Monte Carlo simulation. This can be done using data provided by the cargo insurance companies, which collected the accident occurrences and the cargo losses involved (Table 4).

Figure 5 illustrates an example of how the probabilities are distributed. The maximum cargo value for each range and the cost to the road freight company is considered as being 1% of its value; a deductible that cargo insurances usually charge. Thus, in the occurrence of an accident in the range between $0.01 and R$200,000.00, the value to be considered is always the highest of the range, which is in this case $200,000.00, thus the deductible cost for the carrier should be $2,000.00.

Figure 5 also shows the accumulated percentage values, between 0 and 1, for each accident cost on the right. This is important for the Monte Carlo simulation, which at each iteration selects a random value between 0 and 1 that corresponds to an accident cost. For example, according to Figure 5, each percentage range is equivalent to its cost, and therefore any value selected in the range between 0 and 0.990971 will correspond to an accident cost equal to $0.00.

Thus, a number of 1,000,000 iterations are performed for each arc i j of the problem and the average risk cost of an accident r _{i j} is estimated and it is used in the objective function of the CVRP. This number of iterations was chosen because the accident probabilities are very low, and thus the values of r _{i j} presented a better convergence. However, when the number of iterations is greater, the resolution time also increases, but r _{i j} results show few differences. Therefore, for this problem, the amount of 1,000,000 iterations was ideal.

Figure 5
Accident probabilities and their costs in the Monte Carlo simulation.

Bi-objective in the CVRP

The Capacitated Vehicle Routing Problem (CVRP) was used in a real problem of a freight transportation company located in the city of Limeira/SP. The connection among the cities in the region was represented by graph G = (V, A), where V = 0, ..., n is the set of vertexes representing the cities and A is the set of arcs between i and j. In the problem, there are k identical vehicles each one with capacity cap, and they start and finish at a depot in Limeira. The CVRP minimizing logistics cost (c _{i j} ) and risk cost (r _{i j} ) was also tested in four different instances in which the number of cities varied between 10 to 18 and the number of trucks ranged from 3 to 8, and this is because the number of cities increased, the demand also changed. The instances were named as follows: n10 −k3, n12 −k4, n15 −k8, n18 −k8.

The model is a two-index vehicle flow formulation, as described in Toth & Vigo (2002), where X _{i j} is a binary variable which takes value 1 if a vehicle traverses an arc (i, j) ∈ A and takes value 0 otherwise. The Mixed Integer Programming model was implemented in Python and used Gurobi Optimizer (version 10.0.0) to solve the CVRP.

The PROMETHEE II method

PROMETHEE II is a Multi Criteria Decision Making (MCDM) approach designed to deal with conflicting criteria. The method consists of ranking alternatives through a pairwise comparison by different weighted criteria and it comprises four steps, as discussed in Brans & Smet (2016).

• Step 1: A preference function P _ω (a _a , a _b ) is built to compare alternatives g _ω (a _a ) and g _ω (a _b ) and takes 1 if a _a is preferable to a _b in a criteria ω, (ω ∈ Ω), or 0 otherwise, as Expression (10).

$P_{\omega }(a_a,a_b)=\left\{\begin{array}{l}1,\text{if}{g}_{\omega }\left(a_a\right)<g_{\omega }\left(a_b\right)\\ 0,\text{if}{g}_{\omega }\left(a_a\right)\ge g_{\omega }\left(a_b\right)\end{array}\right.$ (10)

• Step 2: The sum of preference (π) of each criteria ω weighted by w _ω is described by Expression (11)

$\pi (a_a,a_b)=\sum _{\omega \in \Omega }{P}_{\omega }(a_a,a_b)w_{\omega },\forall (a,b)\in M,a\ne b$ (11)

• Step 3: The positive ranking flow (ϕ+) and the negative ranking flow (ϕ ^- ) described by Expressions (12) and (13), respectively, represents how alternative a _a is ranked compared to all others. M is the set all alternatives.

${\varphi }^{+}\left(a_a\right)=\frac{1}{\left|M\right|-1}\sum _{b\in M}\pi (a_a,a_b),\forall a\in M$ (12)

${\varphi }^{-}\left(a_a\right)=\frac{1}{\left|M\right|-1}\sum _{b\in M}\pi (a_b,a_a),\forall a\in M$ (13)

$\varphi \left(a_a\right)={\varphi }^{+}\left(a_a\right)-{\varphi }^{-}\left(a_a\right)$ (14)

• Step 4: Finally, the ranking flow ϕ(a _a ) is a difference of ϕ+(a _a ) and ϕ−(a _a ) as shown in Expression (14) where the higher the ϕ(a _a ), the better the alternative (a _a ) compared to others.

In this paper, a _a and a _b are the arcs of the problem and ω is the criteria related to logistics and risk costs. The weight w _ω ∈ [0, 1] and if w ₁ increases as a proportion of p, then w ₂ decreases (1 − p), once $\sum _{\omega \in \Omega }{w}_{\omega }=1$. As the alternatives are the arcs and they are symmetric, one can change an alternative a _a (∀a ∈ M) to an arc (i, j) (i.e. ϕ(a _a ) = ϕ _ij ), thus ϕ _ij is used as the parameter in the objective function of the CVRP as described below:

$min\sum _{i\in V}\sum _{j\in V}{\varphi }_{ij}.X_{ij}$ (15)

$\sum _{i\in V}{\mathit{\text{X}}}_{ij}=1\forall j\in V\backslash \left\{0\right\}$ (16)

$\sum _{j\in V}{\mathit{\text{X}}}_{ij}=1\forall i\in V\backslash \left\{0\right\}$ (17)

$\sum _{i\in V}{\mathit{\text{X}}}_{i0}=K$ (18)

$\sum _{j\in V}{\mathit{\text{X}}}_{0j}=K$ (19)

${\mathit{\text{U}}}_i-{\mathit{\text{U}}}_j+cap.{\mathit{\text{X}}}_{ij}\le cap-d_j\forall (i,j)\in V\backslash \left\{0\right\}i\ne j|d_i+d_j\le cap$ (20)

$d_i\le {\mathit{\text{U}}}_i\le cap\forall i\in V\backslash \left\{0\right\}$ (21)

${\mathit{\text{X}}}_{ij}\in \{0,1\}\forall (i,j)\in V$ (22)

Expression (15) is the objective function that minimises the ranking flow (ϕ _ij ). The objective function value is accounted for in the variable Z. Expressions (16) and (17) impose that for each vertex representing a location there is only one entry and one exit, respectively. Expressions (18) and (19) represent the depot where it receives and leaves a number of arcs equal to the number of vehicles k. Expressions (20) and (21) impose the capacity limit cap and the connectivity requirements along the route, where Constraints (20) guarantee the sub-tour elimination proposed by Miller et al. (1960) and adapted to CVRP. The continuous variable U _i represents the load of the vehicle after visiting customer i(Toth & Vigo, 2002). The domain of the decision variables of the routes are described in Expression (22).

The ε-constrained method

The ε-constrained method optimises just one objective taking the others as constraints by limiting them to an upper bound defined by parameter ε. In the ε-constrained non-dominated efficient solutions can be produced and other advantages are the linearity preservation of the original problem and only one constraint is added to the model in which the best convenient one can be selected.

In this study, the logistic cost is kept as the objective to be minimized and adds the risk cost as the constraint to be met. Therefore, the difference of this method with the previous one is that the objective function Eq. (23) only minimises c _{i j} .X _{i j} and it is added to the set of the remaining constraints of the CVRP with the new one limited to the risk (Eq. (24)). A variation of ε is made within the range of feasible region to obtain different solutions.

$min\sum _{i\in V}\sum _{j\in V}{c}_{ij}.{\mathit{\text{X}}}_{ij}$ (23)

$\sum _{i\in V}\sum _{j\in V}{r}_{ij}.{\mathit{\text{X}}}_{ij}\le ϵ$ (24)

4 EXPERIMENTAL RESULTS

A freight transportation company in Limeira/SP collaborated in this study supplying real data. The set of all locations in the instance is shown on the map in Figure 6, where the depot in Limeira is identified by the letter ”D”.

Calculating the risk costs

The risk cost (r _ij ) was obtained by the Monte Carlo simulation and a part of the results are shown in Table 5. It was observed that r _ij was obtained in line with what was expected. When comparing the arcs that have the same type of roads such as Piracicaba to Santa Barbara d’Oeste and Limeira to Mogi Mirim, the r _ij for the first arc is greater in the risk cost (r _ij = 1339.24) than the second (r _ij = 150.21), and it was expected because of the higher traffic (6,819 heavy vehicles) on the road h = SP304. On the other hand, the second arc is road h = SP147 where the heavy vehicle flow is only 757 making it a low-risk arc.

Figure 6
Set of locations in all instances of this study.

In the arcs, Mogi Mirim to Rio Claro (r _ij = 464.18) and Araras to Paulinia (r _ij = 859.05), the types of roads are equal for both, but in the first one, the truck goes by a longer distance (l _h = 47.6 km) on a safer road (h = SP147) whereas it is the opposite for the second arc. The same occurs when Limeira to Cosmópolis and Araras to Rio Claro are compared. On the second, the truck runs most of its way (l _h = 21.7 km) on a safer road (h = SP191) making the risk cost of this arc lower (r _ij = 345.60) than the first (r _ij = 1082.12), where the truck travels on two roads with high levels of risk.

When comparing arcs that present similar or nearly a vehicle flow such as Mogi Mirim to Araras and Limeira to Mogi Mirim, it is noted that r _{i j} for the first is greater than the second due to the Mogi Mirim to Araras arc, which is built by a road (h = SP191) single-lane two-way road, according to the nomenclature of CNT (2019), which has the highest death rate among all road types.

It can be observed that SP147 is a very safe road because of its road type, two-lane two-way road with a central safety lane, and the low flow of vehicles. When a vehicle travels long distances on roads such as the SP147, r _ij tends to be low and the arc is safe. On the other hand, roads SP304 and SP330 penalise the r _ij due to their high flow of vehicles.

Other analyses, such as this, were also carried out and the same conclusions were reached. Thus, it can be stated that the analytical approach followed to find r _ij was satisfactory as consistent results were found according to those expected.

The risk cost of using a road, that will be used for all arcs, obtained from the Monte Carlo simulation is presented in Table 5.

Solving the bi-objective function in the CVRP

The optimized results of the CVRP with the bi-objective methods are presented in Figure 7. Analysing the Pareto Frontier of each instance, the ε-constrained method resulted in more varied solutions than Promethee II as setting an upper bound for the constraints could lead to finding weak efficient solutions.

However, it is difficult to find an effective ε interval covering all possible feasible solutions, that is, the smallest value of ε should be found that still makes the solution feasible or a value of ε high enough covering all possible optimal solutions. Therefore, in the case of this study and according to Equation (24), the ε must be greater than the minimum risk value, otherwise the solution is unfeasible, or as high as the maximum risk to cover all possible optimal solutions.

Thus, the PROMETHEE II method can be used as a support to find the range of values of ε. For the results of instance n10-k3 according to Figure 7a and Table 6, when w _ω = 0 the risk assumes the maximum value (r _ij .X _ij = 4616.64). On the other hand, when w _ω = 1 the risk is minimum (r _ij .X _ij = 3463.52). Therefore, in instance n10-k3 when ε < 3463.52, the solution is infeasible, and when ε = 4616.64, it is large enough to cover all optimal solutions, making the ideal interval of ε equal to 3463.52 ≤ ε ≤ 4616.64.

To detail the routes of each solution that compose the Pareto frontier of each instance, see Figures 8, 9, 10 and 11. The cost details of each route of each instance are explained in Tables 6, 7, 8 and 9, respectively.

Table 5
Risk cost (ri j ) obtained from Monte Carlo simulation.

Arcs (i j)	cij($)	rij($)	h	Type of Road	Flow of heavy vehicles	lh(km)
Limeira - Cosmópolis	51.70	1082.12	SP330	Two-lane two-way road with central safety lane	4980	16.5
			SP133	Single-lane two-way road	3373	16.4
Holambra - Cosmópolis	50.76	418.44	SP107	Single-lane two-way road	1790	15.2
			SP332	Two-lane two-way road with central safety lane	1842	10.9
Limeira - Mogi Mirim	117.87	150.21	SP147	Two-lane two-way road with central safety lane	757	54.2
Limeira - Piracicaba	80.19	147.80	SP147	Two-lane two-way road with central safety lane	757	36.9
Limeira - Araras	72.11	971.9	SP330	Two-lane two-way road with central safety lane	4980	31
Limeira - Rio Claro	80.76	974.65	SP330	Two-lane two-way road with central safety lane	4980	12.5
			SP310	Two-lane two-way road with central safety lane	4950	24
Araras - Santa Bárbara d’Oeste	92.56	702.63	SP330	Two-lane two-way road with central safety lane	4980	14.3
			SP348	Two-lane two-way road with central safety lane	2490	38.7
			SP304	Two-lane two-way road with central safety lane	6819	5.9
Santa Bárbara d’Oeste - Sumaré	77.71	536.04	SP304	Two-lane two-way road with central safety lane	6819	4.7
			SP348	Two-lane two-way road with central safety lane	2490	24.9
Araras - Rio Claro	40.38	345.60	SP191	Single-lane two-way road	650	21.7
			SP330	Two-lane two-way road with central safety lane	4980	4
Mogi Mirim - Rio Claro	177.26	464.18	SP147	Two-lane two-way road with central safety lane	535	47.6
			SP330	Two-lane two-way road with central safety lane	4980	6
			SP310	Two-lane two-way road with central safety lane	4150	23.6
Piracicaba - Santa Barbara d’Oeste	44.78	1339.24	SP304	Two-lane two-way road with central safety lane	6819	28.5
Araras - Paulínia	163.00	859.05	SP330	Two-lane two-way road with central safety lane	4980	35.1
			SP147	Two-lane two-way road with central safety lane	535	16.5
			SP332	Two lane two-way road with central safety lane	1603	19.9
Piracicaba - Americana	77.47	1349.1	SP304	Two-lane two-way road with central safety lane	4980	49.3
Mogi Mirim - Araras	110.20	217.33	SP191	Single-lane two-way road	650	48.3
Jundiaí - Socorro	286.54	430.29	SP360-1	Single-lane two-way road	1458	22.0
			SP063	Single-lane two-way road	1267	40.0
			SP008	Single-lane two-way road	1301	47.0
Amparo - Jundiaí	128.06	327.31	SP360-2	Single-lane two-way road	585	40.0
			SP065	Two-lane two-way road with central safety lane	5780	4.4
			SP063	Single-lane two-way road	1267	9.5
			SP360-1	Single-lane two-way road	1458	19
Araras - São João da Boa Vista	222.20	498.75	SP330	Two-lane two-way road with central safety lane	4980	38
			SP225	Single-lane two-way road	880	50
			SP344-1	Single-lane two-way road	862	24
Cosmópolis - Amparo	106.39	610.90	SP332	Two-lane two-way road with central safety lane	1842	10.9
			SP107	Single-lane two-way road	1790	47.1
			SP95-2	Single-lane two-way road	1765	9.7

Figure 7
Pareto frontier from optimization with PROMETHEE II and ε-constrained methods.

4.1 Analysis of Solutions for instance n10-k3

Figure 8 presents the routes that correspond to optimal solutions that compose the Pareto Frontier for instance n10-k3 shown in Figure 7a. The total risk cost, total logistic cost and which roads belong to each solution are shown in Table 6.

Table 6
Routes for each solution of instance n10-k3 presented in Figure 8.

Figure	ε	wω	cij Xij ($)	∆cij Xij	rij Xij ($)	∆rij Xij	Routes
8a	4650	0 to 0.45	729.69	-	4616.64	-	1)D-Holambra-Mogi Mirim-Araras-D
							2)D-Piracicaba-Rio Claro-D
							3)D-Santa Bárbara d’Oeste-Sumaré-Paulínia-Cosmópolis-D
8b	4500	-	778.01	↑ 6.62%	4477.14	↓ 3.02%	1)D-Santa Bárbara d’Oeste-Cosmópolis-Paulínia-Sumaré-D
							2)D-Piracicaba-Rio Claro-D
							3)D-Araras-Holambra-Mogi Mirim-D
8c	4250	-	833.67	↑ 14.25%	4238.76	↓ 8.19%	1)D-Mogi Mirim-Araras-Santa Bárbara d’Oeste-D
							2)D-Piracicaba-Rio Claro-D
							3)D-Cosmópolis-Holambra-Sumaré-Paulínia-D
8d	4100	0.5 to 1.0	834.36	↑ 14.34%	3463.52	↓ 24.98%	1)D-Mogi Mirim-Araras-Santa Bárbara d’Oeste-D
							2)D-Piracicaba-Rio Claro-D
							3)D-Sumaré-Paulínia-Cosmópolis-Holambra-D

Figure 8
Optimal solutions from Pareto-Frontier, for instance n10-k3.

In Figure 8b, the green arc from Limeira to Mogi Mirim has been included in a VRP solution when the risk is relevant enough as compared to the initial solution for instances n10-k3 (Figure 8a). The inclusion of green arcs observed in Figures 8b, 8c and 8d provided a greater reduction in the risk cost than the increase in the logistic costs (Table 6). Road SP147 plays an important role in the risk reduction as it is one of the safest roads in the problem.

Limeira to Cosmópolis (red arc in Figures 8a and 8c) is removed from the solution when ε = 4500 due to the high-risk levels of roads SP330 and SP133. Although its logistic cost is low (c _ij = 51.70), the arc returned as a solution when ε = 4250 and removed again when ε = 4100 because the optimization in this case tends to prioritize safety instead of the logistic costs. The same happens for Limeira to Araras (red arc in Figure 8b) which is removed when ε = 4250 and Araras to Santa Barbara d’Oeste (green one in Figure 8c) is considered due to the difference of their risk cost r _ij .

Road SP147 connects Mogi Mirim to Piracicaba, passing by Limeira. When ε = 4500 (Figure 8b), the model selects this road for the vehicle to travel completely due to its high safety level. On the other hand, when the logistics cost is prioritized at ε = 4650 (Figure 8a), a section of this road, the path from Mogi Mirim to Limeira, is not considered as a solution and the Holambra to Cosmópolis arc that comprises the SP107 and SP133 highways, is added.

The model also selects routes composed by SP348, which connects Araras to Santa Bárbara d’Oeste, from ε = 4250 (Figure 8c) due to its low accident risk. Nonetheless, the SP133 that forms the Limeira to Cosmópolis and Santa Bárbara d’Oeste to Cosmópolis arcs, is excluded when the accident risk should be minimum (Figure 8d), as this road is a single-lane two-way road with high truck traffic.

The SP330 road is long and passes by several cities in the region, such as Limeira, Americana, Araras, Sumaré, and it is considered in the solutions only when the logistic cost is prioritized due to heavy traffic, similarly to SP133, which is shorter in distance and more risky for accidents. On the other hand, the SP147 and SP348 roads are preferred when the low accident risk is recommended.

4.2 Analysis of Solutions for instance n12-k4

Figure 9 presents the routes that correspond to optimal solutions that comprise the Pareto Frontier of instance n12-k4 (Figure 7b). Table 7 summarizes the key results.

The Limeira to Mogi Mirim, Limeira to Piracicaba and Santa Bárbara d’Oeste to Sumaré arcs are introduced as a solution when ε = 9900 (green arcs in Figure 9b) and those are kept as solutions at ε = 7750 and ε = 7300. Limeira to Mogi Mirim and Limeira to Piracicaba are considered because these arcs are formed by the SP147 road, which is a very safe way, as already mentioned in the previous instance. Santa Bárbara d’Oeste to Sumaré are also included as the SP348 road that makes up this arc has low truck traffic and it is a two-lane two-way road with a central safety lane.

Table 7
Routes for each solution of instance n12-k4 presented in Figure 9

Figure	ε	wω	cij Xij ($)	∆cij Xij	rij Xij ($)	∆rij Xij	Routes
9a	10200	0.0 to 0.20	1535.16	-	9979.97		1)D-Santa Bárbara d’Oeste-Piracicaba-Rio Claro-Araras-D
							2)D-Amparo-Paulínia-D
							3) D-Cosmópolis-Mogi Mirim-São João da Boa Vista-D
							4) D-Americana-Sumaré-D
9b	9900	-	1572.56	↑ 2.04%	8449.35	↓ 15.34%	1)D-Santa Bárbara d’Oeste-Sumaré-D
							2)D-Cosmópolis-Amparo-Mogi Mirim-D
							3) D-Piracicaba-Rio Claro-Araras-São João da Boa Vista-D
							4) D-Americana-Paulínia-D
9c	7550	-	1696.06	↑ 8.77%	7728.52	↓ 22.56%	1)D-Santa Bárbara d’Oeste-Sumaré-D
							2)D-Mogi Mirim-Cosmópolis-Amparo-D
							3)D-Piracicaba-Rio Claro-Araras-São João da Boa Vista-D
							4)D-Paulínia-Americana-D
9d	7300	0.55 to 1.0	1835.59	↑ 16.37%	6953.59	↓ 30.32%	1)D-Santa Bárbara D’Oeste-Sumaré-D
							2)D-Mogi Mirim-Cosmópolis-São João da Boa Vista-D
							3)D-Americana-Paulínia-D
							4)D-Piracicaba-Rio Claro-Araras-Amparo-D

Figure 9
Optimal solutions from Pareto-Frontier for instance n12-k4.

The Piracicaba to Santa Bárbara d’Oeste and Limeira to Araras arcs (red arcs in Figure 9a) have a low logistic cost, and thus are selected as a route when ε = 10200. However, when ε = 9900, ε = 7750 and ε = 7300, these arcs are excluded as a solution due to their high risk cost. The Piracicaba to Santa Bárbara d’Oeste and Limeira to Araras arcs are formed by roads SP304 and SP330, respectively, which have the highest accident risk among all the roads in the problem.

The Araras to São João da Boa Vista and Cosmópolis to Amparo arcs (red arcs in Figure 9c) do not have such a high accident risk as the Piracicaba to Santa Bárbara d’Oeste arc, and therefore they are selected as a solution for the routes at ε = 9900 and ε = 7750. However, when the safety level is even more prioritized at ε = 7300, these arcs are eliminated.

Thus, similarly as the previous solution for instance n10-k3, the routes with the SP147 and SP348 roads are preferable and the routes with the SP304 and SP330 roads are excluded when the safety level is prioritized.

4.3 Analysis of Solutions for instance n15-k8

Figure 10 presents the routes that correspond to optimal solutions that compose the Pareto Frontier of instance n15-k8 (Figure 7c). Table 8 summarizes the key results.

Figure 10a displays the optimal solution with the lowest total logistic cost and consequently with a higher accident risk. This solution includes the Americana to Piracicaba arc (red arc in Figure 10a) which has a low logistic cost, but its accident risk is high due to the high traffic on road SP304. When the accident risk is reduced, other routes are configured as shown in Figures 10b and 10c, where the safe roads are introduced (SP147 and SP191, respectively) as their logistical costs are not high.

The Limeira to Rio Claro arc (Figure 10b) also has low logistics cost and high accident risk because it contains two roads that are dangerous shown by the data, SP310 and SP330. However, this arc is still considered as a solution when the accident risk starts to be prioritized at ε = 13300. However, when the safety level is even more prioritized at ε = 12150 and ε = 12000, the model excludes this arc as a solution and selects a safer one, such as the arc passing by Araras to Rio Claro (Figure 10c).

Figure 10
Optimal solutions from Pareto-Frontier for instance n15-k8.

Table 8
Routes for each solution of instance n15-k8 presented in Figure 10.

Figure	ε	wω	cij Xij ($)	∆cij Xij	rij Xij ($)	∆rij Xij	Routes
10a	15400	0.0	3081.24	-	15374.34	-	1)D-Mogi Mirim-D
							2)D-Cosmópolis-D
							3)D-Santa Bárbara d’Oeste-Tambaú-D
							4)D-Sumaré-Amparo-D
							5)D-Jundiaí-Araras-D
							6)D-São João da Boa Vista-D
							7)D-Paulínia-Pedreira-D
							8)D-Americana-Piracicaba-Rio Claro-D
10b	13300	-	3123.56	↑ 1.37%	13255.95	↓ 13.78%	1) D-Tambaú-D
							2)D-Mogi Mirim-D
							3)D-Amparo-Sumaré-Santa Bárbara d’Oeste-D
							4)D-Cosmópolis-D
							5)D-Jundiaí-Araras-D
							6)D-São João da Boa Vista-D
							7)D-Paulínia-Pedreira-D
							8)D-Americana-Rio Claro-Piracicaba-D
10c	12150	0.05 to 0.40	3260.07	↑ 5.80%	12113.77	↓ 21.21%	1)D-Tambaú-D
							2)D-Mogi Mirim-D
							3)D-Amparo-D
							4)D-Americana-Jundiaí-D
							5)D-Santa Bárbara d’Oeste-Sumaré-Cosmópolis-D
							6)D-São João da Boa Vista-D
							7)D-Paulínia-Pedreira-D
							8)D-Araras-Rio Claro-Piracicaba-D
10d	12000	0.45 to 1.0	3370.42	↑ 9.39%	11758.74	↓ 23.52%	1)D-Tambaú-D
							2)D-Mogi Mirim-D
							3)D-Amparo-D
							4)D-Americana-Paulínia-D
							5)D-Santa Bárbara d’Oeste-Sumaré-Cosmópolis-D
							6)D-São João da Boa Vista-D
							7)D-Jundiaí-Pedreira-D
							8)D-Araras-Rio Claro-Piracicaba-D

4.4 Analysis of Solutions for instance n18-k8

Figure 11 presents the routes that correspond to optimal solutions that compose the Pareto Frontier of instance n18-k8 (Figure 7d). Table 9 summarizes the key results.

The Limeira to Rio Claro arc (red arc in Figure 11a) is introduced as an optimal solution when ε = 14850 due to its low logistics cost. However, when risk is prioritized in ε = 12550, ε = 11700 and ε = 11100, this arc is excluded due to the high accident risk in SP330 and SP310, similarly as indicated in instance n15-k8. The Araras to Rio Claro arc (green arc in Figure 11b) is preferred as a solution when ε = 12550, ε = 11700 and ε = 11100 as it has low risks of accident.

The Jundiaí to Socorro arc (red arc in Figure 11a) is selected as an optimal solution only in ε = 14850. When the risk is prioritized, the routes obtained in ε = 12550, ε = 11700 and ε = 11100 do not contain the Jundiaí to Socorro arc and the new solution introduces the Amparo to Jundiaí arc (green arc in Figure 11b) as option. This is because on the Amparo to Jundiaí arc, the vehicle travels mostly on the SP360-2 road, where the volume of vehicles is low, while on the Jundiaí to Socorro arc, the vehicle travels only on roads with high traffic. Thus, the risk cost for the Amparo to Jundiaí arc is lower and preferable as an optimal solution.

Figure 11
Optimized routes for instance n18-k8.

5 FINAL REMARKS

The need to consider accident risk in VRP is important because when an accident occurs not only are people exposed to risk, but the loss of transported goods can have major consequences such as: financial, interruptions in the supply chain and social-environmental impacts.

This study considered the issue of route safety in the VRP for a cargo transportation company to support decision-makers to choose the best routes that minimise both the logistics costs and accident risks.

Due to the limited data availability of road accidents, an analytical approach based on simple statistic calculations was developed to estimate the accident probabilities for each arc between locations, and the costs related to the accident were estimated through Monte Carlo simulations. The PROMETHEE II and ε-constrained methods were implemented to deal with the bi-objective (i.e. tradeoff between logistic cost and accident risk cost) in the VRP. The obtained results were coherent for the analytical approach as a whole, converging to what was expected of the problem.

Table 9
Routing solution of instance n18-k8 presented in Figure 11.

Figure	ε	wω	cij Xij ($)	∆cij Xij	rij Xij ($)	∆rij Xij	Routes
11a	14850	-	2967.41	-	14836.43	-	1)D-Mogi Mirim-Araras-D
							2)D-São João da Boa Vista-Tambaú-D
							3)D-Cosmópolis-D
							4)D-Santa Bárbara d’Oeste-Sumaré-Americana-Limeira
							5)D-Piracicaba-São Pedro-Rio Claro-D
							6)D-Paulínia-D
							7)D-Pedreira-Amparo-Serra Negra-D
							8)D-Jundiaí-Socorro-D
11b	12550	0.35 to 0.40	3221.46	↑ 8.56%	12544.98	↓ 15.44%	1)D-Mogi Mirim-Serra Negra-D
							2)D-Tambaú-São João da Boa Vista-D
							3)D-Jundiaí-Amparo-Socorro-D
							4)D-Santa Bárbara d’Oeste-Sumaré-D
							5)D-Piracicaba-D
							6)D-Paulínia-D
							7)D-São Pedro-Rio Claro-Araras-Pedreira-D
							8)D-Cosmópolis-Americana-D
11c	11700	-	3542.54	↑ 19.38%	12113.77	↓ 21.46%	1)D-Mogi Mirim-Serra Negra-D
							2)D-Tambaú-São João da Boa Vista-D
							3)D-Jundiaí-Amparo-D
							4)D-Santa Bárbara d’Oeste-Sumaré-D
							5)D-Piracicaba-D
							6)D-Paulínia-D
							7)D-Pedreira-Araras-Rio Claro-São Pedro-D
							8)D-Socorro-Cosmópolis-Americana-D
11d	11100	0.80	3951.87	↑ 33.18%	11090.86	↓ 25.25%	1)D-Mogi Mirim-Serra Negra-D
							2)D-Americana-Cosmópolis-Socorro-D
							3)D-Tambaú-Paulínia-D
							4)D-Amparo-Jundiaí-D
							5)D-São João da Boa Vista-D
							6)D-São Pedro-Rio Claro-Araras-Pedreira-D
							7)D-Santa Bárbara d’Oeste-Sumaré-D
							8)D-Piracicaba-D

Analysing the results of instance n10-k3 presented in Table 6, the increase in logistic cost is greater than the risk cost reduction at ε = 4500, but when the security level of the routes receives higher priority, the risk cost reduction is greater at ε = 4100. On the other hand, in instance n18-k8 (results in Table 9), the opposite happened. In the first moment, the risk cost reduction is greater at ε = 12550 and ε = 11700, but the increase in logistics cost is greater when security is on the highest prioritization level (ε = 11100). For instances n12-k4 (Table 7) and n15-k8 (Table 8), the risk cost reduction is greater than the increases in logistics costs for all safety levels.

Thus, a decision-maker could analyse the total logistic cost and risk of accidents, and one could plan the routes based on which safety level it wishes to operate, as the level of operational safety depends on some factors such as the cargo type and the cargo value.

It can be affirmed that the analytical approach worked coherently by representing the trade-off between logistics costs and the risk of accident due to more dangerous roads. According to the parameters that regulate the bi-objective approaches, the model may prioritize safer routes or the logistic costs, showing different solutions to the decision-maker.

The limitation of the study was the small amount of data and their types to estimate the risk cost in an even better way, as more data and different types of information on road hazards could be used to develop better risk models. For further studies, the analytical approach will be extended to include more conflicting objectives to the logistical cost, such as environmental costs, making the VRP more multi-objective.

References

• ANDROUTSOPOULOS KN & ZOGRAFOS KG. 2012. A bi-objective time-dependent vehicle routing and scheduling problem for hazardous materials distribution. EURO Journal on Transportation and Logistics, 1(1-2): 157-183.
• BILATO GA. 2022. GabrielBilato/paper riskVRP 2022. Available at: https://github.com/GabrielBilato/paper riskVRP2022
  » https://github.com/GabrielBilato/paper riskVRP2022
• BRAEKERS K, RAMAEKERS K & VAN NIEUWENHUYSE I. 2016. The vehicle routing problem: State of the art classification and review. Computers and Industrial Engineering, 99: 300-313.
• BRANKO MILOVANOVIĆ. 2012. Methodology for establishing the routes for transportation of dangerous goods on the basis of the risk level Case study: City of Belgrade. Scientific Research and Essays, 7(1).
• BRANS JP & SMET YD. 2016. PROMETHEE methods. International Series in Operations Research and Management Science, 233: 187-219.
• BULA GA, GONZALEZ FA, PRODHON C, AFSAR HM & VELASCO NM. 2016. Mixed Integer Linear Programming Model for Vehicle Routing Problem for Hazardous Materials Transportation. IFAC-PapersOnLine, 49(12): 538-543.
• BULA GA, MURAT AFSAR H, GONZÁLEZ FA, PRODHON C & VELASCO N. 2019. Bi-objective vehicle routing problem for hazardous materials transportation. Journal of Cleaner Production, 206: 976-986.
• CARRESE S, CUNEO V, NIGRO M, PIZZUTI R, ARDITO CF & MARSEGLIA G. 2022. Optimization of downstream fuel logistics based on road infrastructure conditions and exposure to accident events. Transport Policy, 124: 96-105. Available at: https://www.sciencedirect.com/ science/article/pii/S0967070X19305025
  » https://www.sciencedirect.com/ science/article/pii/S0967070X19305025
• CHAI H, HE R, KANG R, JIA X & DAI C. 2023. Solving Bi-Objective Vehicle Routing Problems with Driving Risk Consideration for Hazardous Materials Transportation. Sustainability, 15(9). Available at: https://www.mdpi.com/2071-1050/15/9/7619
  » https://www.mdpi.com/2071-1050/15/9/7619
• CNT. 2019. Acidentes Rodoviários: Estatísticas envolvendo caminhões. Confederação Nacional de Transportes, a. lccessed: 2021-11-15.
• DER-SP. 2021. Contagem Volumétrica Classificatória VDM. A. lccessed: 2021-11-15. DNIT. 2021. Plano Nacional de Contagem de Tráfego. A. lccessed: 2021-11-15.
• DU J, LI X, YU L, DAN R & ZHOU J. 2017. Multi-depot vehicle routing problem for hazardous materials transportation: A fuzzy bilevel programming. Information sciences, 399: 201-218.
• ENGSTRÖM R. 2016. The Roads’ Role in the Freight Transport System. Transportation Research Procedia, 14: 1443-1452.
• ERKUT E & INGOLFSSON A. 2005. Transport risk models for hazardous materials: Revisited. Operations Research Letters, 33(1): 81-89.
• GHANNADPOUR SF & ZANDIYEH F. 2020. An adapted multi-objective genetic algorithm for solving the cash in transit vehicle routing problem with vulnerability estimation for risk quantification. Engineering Applications of Artificial Intelligence, 96(September): 103964.
• HOLECZEK N. 2021. Analysis of different risk models for the hazardous materials vehicle routing problem in urban areas. Cleaner Environmental Systems, 2(February): 100022.
• MILLER CE, TUCKER AW & ZEMLIN RA. 1960. Integer Programming Formulation of Traveling Salesman Problems. J. ACM, 7(4): 326-329. Available at: https://dblp.org/rec/journals/jacm/ MillerTZ60.bib
  » https://dblp.org/rec/journals/jacm/ MillerTZ60.bib
• PRADHANANGA R, TANIGUCHI E, YAMADA T & QURESHI AG. 2014. Bi-objective decision support system for routing and scheduling of hazardous materials. Socio-Economic Planning Sciences, 48(2): 135-148.
• QUALP. 2021. A. lccessed: 2021-10-12.
• TALARICO L, SÖ RENSEN K & SPRINGAEL J. 2015. Metaheuristics for the risk-constrained cash-in-transit vehicle routing problem. European Journal of Operational Research, 244(2): 457-470.
• TALARICO L, SPRINGAEL J, SÖ RENSEN K & TALARICO F. 2017. A large neighbourhood metaheuristic for the risk-constrained cash-in-transit vehicle routing problem. Computers and Operations Research, 78: 547-556.
• TOTH P & VIGO D. 2002. An Overview of Vehicle Routing Problems. In: TOTH P & VIGO D (Eds.), The Vehicle Routing Problem, vol. 9 of SIAM monographs on discrete mathematics and applications. pp. 1-26. SIAM. Available at: https://dblp.org/rec/books/siam/02/TothV02.bib
  » https://dblp.org/rec/books/siam/02/TothV02.bib
• WANG N, ZHANG M, CHE A & JIANG B. 2018. Bi-Objective Vehicle Routing for Hazardous Materials Transportation with No Vehicles Travelling in Echelon. IEEE Transactions on Intelligent Transportation Systems, 19(6): 1867-1879.
• ZHENG B. 2010. Multi-objective vehicle routing problem in hazardous material transportation. ICLEM 2010: Logistics for Sustained Economic Development Infrastructure, Information, Integration Proceedings of the 2010 International Conference of Logistics Engineering and Management, 387: 3139-3145.

Publication Dates

• Publication in this collection
  10 Nov 2023
• Date of issue
  2023

History

• Received
  13 Dec 2022
• Accepted
  30 July 2023