Fuzzy green vehicle routing problem for designing a three echelons supply chain

Journal Pre-proof Fuzzy green vehicle routing problem for designing a three echelons supply chain

Antonio Giallanza, Gabriella Li Puma PII:

S0959-6526(20)30821-0

DOI:

https://doi.org/10.1016/j.jclepro.2020.120774

Reference:

JCLP 120774

To appear in:

Journal of Cleaner Production

Received Date:

12 June 2019

Accepted Date:

25 February 2020

Please cite this article as: Antonio Giallanza, Gabriella Li Puma, Fuzzy green vehicle routing problem for designing a three echelons supply chain, Journal of Cleaner Production (2020), https://doi.org/10.1016/j.jclepro.2020.120774

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier.

Journal Pre-proof

Fuzzy green vehicle routing problem for designing a three echelons supply chain Antonio Giallanza*, Gabriella Li Puma Engineering Department, University of Palermo, Viale delle Scienze, Palermo, Italy *Corresponding author: [email protected]

Journal Pre-proof Abstract In this study, a three-echelon fuzzy green vehicle routing problem (3E-FGVRP) is considered for designing a regional agri-food supply chain on a time horizon. To account for the variability associated with the quantities requested by customers, it is assumed that the demands are fuzzy numbers simulated by a time-dependent algorithm. Moreover, the vehicle fleet and distribution centres are considered with a defined capacity. The credibility theory of fuzzy sets is used to implement a multi-objective fuzzy chance-constrained programming model, where the total costs and carbon emissions are minimised. The resolution of the 3E-FGVRP is conducted by using a nondominated sorting genetic algorithm. The multiple-criteria decision-making ELECTRE III method is applied to select the best solutions belonging to each Pareto front. Finally, the validity of the model is demonstrated by performing an optimisation procedure with three different initial random sets of populations. The application of the model to a case study of the Sicilian agri-food context confirms the robustness of the model, and the optimal configurations of the three-echelon supply chain can be found. Keywords: GVRP, simulation, fuzzy demand, credibility theory, Multi objectives optimization, NSGA-II 1. Introduction and literature review The design of a distribution network is an important issue that needs to be addressed because it influences the success of companies. The emergence of new business models that are closer to the needs and lifestyles of consumers has led to the development of new industrial and, for the most part, distributive realities. The huge investments in marketing and the costs necessary for the distribution, as well as for the acquisition of a fleet of vehicles, make it necessary to minimise the costs of transportation and, more generally, to manage routes. If companies do not correctly manage these aspects, they would be forced to increase the delivery costs or in the worst case scenario, the price of the goods, to achieve higher revenues. Both scenarios are very risky because they would increase the possibility of customers purchasing substitute products, thereby decreasing sales volumes and losing customers’ trust. Therefore, it is clear that designing a robust supply chain is essential for a company to survive in the increasingly competitive market. The decisions about the locations of distribution facilities, amount of inbound and outbound flows, and definition of routes represent only some of the innumerable factors to be necessarily considered. To give concrete solutions to these issues, the vehicle routing problem (VRP) has been implemented

Journal Pre-proof for many decades. The VRP aims at determining the optimal routes for a fleet of vehicles under several constraints involving depot capacity, temporal priorities, and variable customer demand. The body of knowledge has a wide literature in the field of VRP, which is briefly described as follows. The first application of VRP was attributed to Dantzing and Ramser [1], who proposed a study concerning the route definition of a fleet of gasoline delivery trucks between a single terminal and a large number of service stations. The paths characterised by the minimum travelled distance are computed by applying a linear programming formulation. After this work, several applications and variants of the VRP have been implemented. References [2–4] are some of the first works dealing with the VRP in the contemporary approach. Orloff [2] defined the general routing problem starting from the travelling salesman problem and the Chinese postman problem, as a peculiar limiting case. The paper proposed a resolution algorithm from which the approach managed to reduce the effective problem size. Christofides [3] dealt with the problem of designing routes for a vehicle fleet by minimising the cost of distribution. The study discussed exact and approximate methods for solving VRPs. Golden et al. [4] applied heuristic algorithms for a VRP by proposing modifications and extensions that allow solving the problem with hundreds of demand points in a few seconds. The increasing transportation cost is a serious concern for many business organisations faced with distribution issues. As a matter of fact, the growing complexity of decision making in the real industrial environment, due to the design of increasingly complex business models, causes more difficulties in the VRP research. Some of these factors include traffic conditions, government regulations, time restrictions, and sustainability aspects. As examples, the following studies concerning temporal issues are presented. Malandraki et al. [5] pioneered a study in which the now famous VRP model is integrated with temporal variables. That study is an important work that allows including time-based features to VRP models. Afterwards, several studies have dealt with these topics. Solomon [6] analysed algorithms for vehicle routing and scheduling problems with timewindow constraints. The problem set includes routing and scheduling environments that differ in terms of the type of data used to generate the problems. Savelsbergh [7] proposed a routing problem with time windows. The local search algorithm with the addition of time windows presents the checking of feasibility constraints that requires an additional time effort. Vidal et al. [8] presented a hybrid genetic search with advanced diversity control for many time-constrained VRPs, introducing several new features to manage the temporal dimension. In particular, new techniques allowing to penalise infeasible solutions and reformulate routes with respect to time-based constraints are described. Another extension of the VRP model, widely applied in the literature, considers input demand data affected by uncertainty. Teodorović et al. [9] proposed a model to solve a VRP by applying a heuristic

Journal Pre-proof sweeping algorithm with the rules of fuzzy arithmetic and logic, when demands at the nodes are uncertain. They considered the routes between one depot from which vehicles depart and return to several customers. Zheng et al. [10] studied a VRP with vehicle travel times as fuzzy variables. The implemented fuzzy optimisation model is a hybrid intelligent algorithm that takes a credibility measure into consideration by integrating a fuzzy simulation and genetic algorithm. Erbao et al. [11] presented a VRP with fuzzy demands (VRPFD) by applying a fuzzy chance-constrained programming model based on the fuzzy credibility theory. They focused on the dispatcher preference index’s impact on the objective value by means of a stochastic simulation. The fuzzy chance-constrained programming model was also used by Cao et al. [12] to solve an open VRPFD. This type of problem differs from the classic one in a key aspect, i.e. the vehicle is not required to return to the distribution depot after servicing the last customer on its route. Brito et al. [13] proposed a close–open VRP, the nearest model to the outsourcing distribution, which is a very usual practice in which delivery companies are paid based on the distance covered by their vehicles. A hybrid metaheuristic algorithm is applied by considering the customer demands and travel times at imprecise data; thus, the capacity and time-window constraints are considered as fuzzy constraints. Shi et al. [14] studied a vehicle routing scheduling problem related to home health care (HHC) companies in which the optimal routes are found by minimising the distribution total costs by using a hybrid genetic algorithm integrated with stochastic simulation methods. In particular, the amount of delivered drugs to each patient is considered non-deterministic, according to a survey of the HHC companies; thus, the programming model is developed by considering drug demand as a fuzzy variable. Environmental disasters afflicting the planet, caused mainly by greenhouse gas emissions, have prompted international and European legislators to impose restrictions in this context, causing companies to limit emissions and invest in renewable energies. For this reason, scientific research related to logistics fields has also moved in this direction, developing a new branch of models, namely green vehicle routing problem (GVRP). The above-mentioned models of VRP are traditionally designed to minimise the delivery process costs in terms of vehicle fleet size or their travelled distances. The sustainable environmental extension of the problem routing has led to the introduction of other objective criteria for route determination, often added to the cost function, by which the environmental impact must be decreased. The reasons why these parameters are included in the VRP research do not only reside in the legislative imposition, but also constitute great reputational benefits for companies. Several studies [15–25] showed that environmental issues play an increasingly important role in corporate social responsibility policies and, despite the huge investments required, the acquisition of a green credential has been giving companies a great competitive advantage.

Journal Pre-proof Different applications can be found in the field of GVRPs. Lin et al. [26] provided a classification of the GVRP by identifying three macro-categories: GVRP, pollution routing problem (PRP), and VRP in reverse logistics. Slight differences exist between these problems. The GVRP aims at optimising the energy consumption of transportation. PRPs aim at determining delivery schemes characterised by less pollution (i.e. reduction of carbon emissions). The VRP in reverse logistics focuses on the distribution aspects of reverse logistics, defined in [26] as ‘the process of planning, implementing, and controlling backward flows of raw materials in process inventory, packaging, and finished goods, from a manufacturing, distribution, or use point to a point of recovery or point of proper disposal. Based on the above classification, some of the most significant studies are discussed below. Jemai et al. [27] presented a bi-objective GVRP in which the routes for serving a set of customers are found by minimising the total travelled distance of vehicles and the CO2 emission. The problem is solved by applying the evolutionary non-dominated sorting genetic algorithm (NSGA-II) and performing benchmark and statistical analyses to validate the results. The same objective functions are considered in [28]; however, the authors used a block recombination approach as the optimisation method. Different clusters for each city visited by different trucks are generated and each of these is approximated to a block. The implemented method is then compared to a conventional evolutionary algorithm to compare the experimental results. Koç et al. [29] and Erdoĝan et al. [30] studied a GVRP by considering a limited number of alternative fuel-powered vehicles and limited refuelling infrastructure. The resolution model in [29] is developed by means of a simulated annealing heuristic based on an exact solution approach, whereas in [30], a mixed-integer linear program and two heuristics (i.e. the modified Clarke and Wright savings heuristic and the density-based clustering algorithm) are developed. Bektaş et al. [31] presented a PRP by taking into consideration the amount of greenhouse emission, fuel, travel time and cost, as well as travel distance. They implemented models with and without time windows and finally validated them by applying realistic instances. Demir et al. [32] applied an adaptive large neighbourhood search for a PRP, which involves determining the speed of vehicles that serve a set of customers on each route segment to minimise a function comprising fuel, emission, and driver costs. Kramer et al. [33] proposed a single-objective model, where operational and environmental costs (the latter based on driver’s wages and fuel consumption) are minimised. The proposed optimisation pattern combines a local search-based metaheuristic with an integer programming approach that enables a tighter integration of route and speed decisions. Finally, Hosseini-Nasab et al. [34] dealt with a PRP by considering three objective functions, two of which have to be minimised (trip cost and fuel consumption) and the third has to be maximised, and involves

Journal Pre-proof addressing customer satisfaction for delivering goods in a shorter time. They considered a penalty of fuel consumption for deciding which feasible arches have to belong to routes. The VRP in reverse logistics category has a wider literature; more themes can be addressed to it according to Lin et al. [26]: selective pickups with pricing, waste collection, end-of-life goods collection, simultaneous distribution, and collection. In Micale et al. [35] in order to satisfy the complicated limits of the real world different types of vehicles in terms of maximum capacity, velocity and emissions, asymmetric paths, vehicle-client constraints and delivery time windows were considered. In particular the authors propose a firefly algorithm to resolve the VRP problem and to consider simultaneously different aspects (e.g. economic and environmental) use the TOPSIS technique. The category of selective pickup and delivery problem addresses real-world issues in the logistics industry and aims at finding the shortest route that can supply delivery nodes with required commodities from some pickup nodes. This problem can substantially reduce the transportation cost and fits real-world logistics scenarios. Ting and Liao [36] proposed a memetic algorithm based on a genetic algorithm with local search to resolve the selective pickup and delivery problem by proving, furthermore, that it is an NP-hard problem. They showed a novel representation of candidate solutions for the selection of pickup nodes. Others studies concerning this topic are reported in [37–45]. One of the key areas of reverse logistics concerns waste management. The majority of manufacturing processes involve generation of waste materials and products that cannot be often disposed of inhouse. Thus, the need for organising and managing the collection of such materials in a cost-effective and ecological manner arises. The waste collection concerns are tackled in [46–53]. For example, Hannan et al. [46] studied a modified particle swarm optimisation algorithm in a capacitated VRP model to determine the best waste collection and route optimisation solutions. The study proved that the more the waste collection route is optimised, the more reduction in different costs and environmental effects there will be. Several studies have analysed the end-of-life goods collection area [54–58]. Among them, Habibi et al. [57] showed a joint optimisation of collection and disassembly processes to improve the global performance of a reverse supply chain of waste electrical and electronic equipment. They presented a collection of end-of-life products dealing with a VRP and solved it by including a lower total cost corresponding to the component demand satisfaction. Finally, some studies [59–61] solved the routing problems of reverse logistics in the field of simultaneous distribution and collection. Tasan and Gen [61] applied a genetic algorithm-based approach to solve a VRP with simultaneous pickups and deliveries. They considered a fleet of vehicles originating from a depot that serves customers with pickups and deliveries from/to their locations.

Journal Pre-proof As mentioned above, each study applies one or more optimisation algorithms. For the sake of completeness, other relevant works not mentioned are classified according to the implemented optimisation. In particular, references [62–67] belong to the body of knowledge in which VRPs are solved by means of exact optimisation algorithms. On the other hand, references [68–79] dealt with the application of heuristic and metaheuristic optimisation methods, among which, the genetic algorithm, simulated annealing, and particle swarm optimisation are the most commonly used. The literature review shows that heuristic approaches are preferred in most recent studies rather than exact methods. This study shows the application of NSGA-II to solve a fuzzy GVRP for a three-echelon (3E-FGVRP) supply chain under uncertain customer demand. The mathematical programming model, which is completely developed in the MATLAB environment, takes into consideration triangular fuzzy numbers as input data; thus, the credibility theory is used to solve the fuzzy chance-constrained programming model. Moreover, the 3E-FGVRP aims at proving the feasibility in terms of the costs involved and environmental impact of a delivery supply chain. This study serves as a preliminary analysis for the implementation of a delivery agri-food supply chain that is aimed at enhancing the local small and medium-sized enterprises in the agri-food sector and, at the same time, allowing inhabitants to consume real zero-mile regional agricultural products. Being an introductory study, the demand values related to each customer group are simulated by developing a customised programming code in which the temporal variations in customers’ purchasing behaviour are taken into consideration. To the best of the authors’ knowledge, the demand simulation algorithm has not yet been considered in VRPs; thus, this study proposes an innovative tool to ascertain the feasibility of a distribution chain under the conditions of uncertainty. The mathematical programming model includes two objective functions to be minimised: the total costs involved in the supply chain and the CO2 emissions of the vehicle fleet. The resolution of the aforementioned model is conducted by means of a heuristic NSGA-II related to the 3E-FGVRP by which the Pareto fronts related to each time period are found. Therefore, to choose the respective best solutions, the multiple-criteria decisionmaking (MCDM) ELECTRE III is applied. Finally, a sensitivity analysis is conducted to test the robustness of the implemented optimisation algorithm with different initial random populations. Figure 1 shows the main steps of the implemented methodology:

Journal Pre-proof

Triangular fuzzy data

ELECTRE III methodology

NSGA

a. criteria importance

INPUT DATA

GOAL GVRP fuzzy

method

PARETO FRONT

SENSITIVITY ANALYSIS

b. thresholds definition

c. Final ranking

Figure 1: main steps of the implemented methodology To facilitate reading this paper, its structure is presented as follows. Section 2 briefly describes the most significant and useful notions of the fuzzy credibility theory. The formulation of the fuzzy chance-constrained programming model, customer demand simulation, and the proposed NSGA-II are discussed in Sections 3, 4, and 5, respectively. The ELECTRE III method is described in Section 6 and finally, the computational results and sensitivity analysis are presented in Section 7. 2. Fuzzy credibility measure theory The concept of fuzzy sets was first addressed by Zadeh [80] and, since then, it has been applied to a wide variety of real problems. The fuzzy theory has become a useful tool for managing problems that are strongly characterised by uncertainty. Conventionally, to study the fuzzy phenomena’s behaviour, the possibility measure theory, also proposed by Zadeh [81], was widely used in the literature and it can be considered as the counterpart of the probability theory for crisp numbers. However, the credibility theory defined by Liu [82] can also be used for dealing with fuzziness in decision systems. To improve the understanding of the proposed methodology, some basic concepts of fuzzy theory and axioms of the possibility measure theory are provided below. Let X be a space of objects, where its generic element is denoted by x, so that X = {x}. A fuzzy set A in X is characterised by a membership function μA(x), which associates a real number in the [0, 1] interval with each point of X, that means the ‘grade of membership’ of x in A. Let Θ be a non-empty set, and let P(Θ) be the power of Θ, where each element in P(Θ) is represented by an event. A nonnegative number Pos{A} is associated with each event A ∈ P(Θ), which represents the possibility that event A will occur. Because ϕ is an empty set, the possibility measures have the mathematical properties indicated in the following axioms [82]:

Journal Pre-proof  Axiom 1. Pos {ϕ} = 0.  Axiom 2. Pos {Θ} = 1.  Axiom 3. ∀ 𝐴𝑘 ∈ P (Θ), Pos{∪k Ak} = supk Pos(Ak), where Ak is an arbitrary collection in P(Θ).  Axiom 4. If Θk is a non-empty set and the set function Posk {}, k = 1, 2, . . . , n, satisfies the conditions stated in the previous axioms, and Θ = Θ 1 × Θ2 × … × Θn, then for each A in P(Θ), Pos{A} = sup(θ1,θ2,...,θn)∈A Pos1{Θ1} ∧ Pos2{Θ2} ∧ … ∧ Posn{Θn}. The credibility measure theory is completely derived from the above four axioms and is summarised in the definitions below [82]. Definition 1. Let (Θ, P(Θ), Pos) and A be a possibility space and a generic set in P(Θ), respectively. The necessity measure of A is defined as Nec {A} = 1 – Pos {A}. Definition 2. Let (Θ, P(Θ), Pos) and A be a possibility space and a generic set in P(Θ), respectively. 1

The credibility measure related to A is defined by the relation Cr {A} = 2 (Pos {A} + Nec {A}). Let us consider a general triangular fuzzy number ã = (a1, a2, a3), where a1 < a2 < a3. The µã(x) expression denotes the membership function whose expression is given by Eq. (1), where a2 is the corresponding value of the highest grade of membership.

𝜇𝑎(𝑥) =

{

𝑥 ― 𝑎1

𝑎2 ― 𝑎1, 𝑎3 ― 𝑥

𝑎1 ≤ 𝑥 ≤ 𝑎2

𝑎3 ― 𝑎2,

(1)

𝑎2 ≤ 𝑥 ≤ 𝑎3 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

0,

The possibility, necessity, and credibility measures indicated in Eqs. (2)–(4) are derived by means of the above-mentioned rules and assumptions [82]. Pos{𝑎 ≥ 𝑏} = 𝑠𝑢𝑝{𝜇𝑎(𝑥)∣𝑥 ∈ 𝑅,𝑥 ≥ 𝑏} 𝑃𝑜𝑠{𝑎 ≥ 𝑏} =

{

1 𝑎3 ― 𝑏 𝑎3 ― 𝑎2

0

𝑖𝑓 𝑏 ≤ 𝑎2 𝑖𝑓 𝑎2 ≤ 𝑏 ≤ 𝑎3 𝑖𝑓 𝑏 ≥ 𝑎3

(2)

Nec{𝑎 ≥ 𝑏} = 1 ― 𝑠𝑢𝑝{𝜇𝑎(𝑥)∣𝑥 ∈ 𝑅,𝑥 < 𝑏} 𝑁𝑒𝑐{𝑎 ≥ 𝑏} =

{

1

1 𝑎2 ― 𝑏 𝑎2 ― 𝑎1

0

𝑖𝑓 𝑏 ≤ 𝑎1 𝑖𝑓 𝑎1 ≤ 𝑏 ≤ 𝑎2 𝑖𝑓 𝑏 ≥ 𝑎2

Cr{𝑎 ≥ 𝑏} = 2[𝑃𝑜𝑠{𝑎 ≥ 𝑏} + 𝑁𝑒𝑐{𝑎 ≥ 𝑏}]

(3)

Journal Pre-proof

𝐶𝑟{𝑎 ≥ 𝑏} =

{

1 2𝑎2 ― 𝑎1 ― 𝑏 2(𝑎2 ― 𝑎1) 𝑎3 ― 𝑏 2(𝑎3 ― 𝑎2)

0

𝑖𝑓 𝑏 ≤ 𝑎1 𝑖𝑓 𝑎1 ≤ 𝑏 ≤ 𝑎2 𝑖𝑓 𝑎2 ≤ 𝑏 ≤ 𝑎3 𝑖𝑓 𝑏 ≥ 𝑎3

(4)

Possibility and necessity measures are not self-dual. This condition implies that an event could occur even if its necessity measure is 0 or it could fail if its possibility assumes a value of 1. To overcome this pitfall, Liu and Liu [83] proposed the credibility fuzzy measure (4) computed as the arithmetic average of the possibility and necessity measures and is characterised by self-duality. The latter property establishes a strong relation between the measure of an event and the measure of its opposite event (each one determines the other one) so that a fuzzy event certainly occurs if its credibility value is equal to 1, and therefore, the event fails if its credibility is 0. 3. Fuzzy chance-constrained programming model for 3E-FGVRP This study presents a preliminary investigation for designing a delivery supply chain in a delimited geographical area. The idea arises from the desire to create a reality that can allow inhabitants to buy and consume zero-mile agricultural products by taking inspiration from the most modern concepts of agri-food distribution. For this purpose, the authors decide to identify nine customer groups in every regional province and to choose food supplier locations that are representative of seven real Sicilian food suppliers classified by product category. On the other hand, the distribution centre locations are chosen in strategic sites of the island for road transportation. The mathematical model that must represent the 3E-FGVRP concept reflects the real dynamics that the supply chain should have once implemented, although considering customer demand as fuzzy data, as this study aims at performing a preliminary design analysis. Furthermore, the authors decide to prove the supply chain in several time periods by considering a quantitative variation in the customer demands that takes into consideration the customers’ real purchasing habits. Figure 2 shows the geographical representation of the supply chain actors, where the squares symbolise food suppliers, the circles represent distribution centres, whereas the triangles indicate the customer groups. In Table 1, the notations of the mathematical model are presented. The distances between each pair of locations are determined through using Google Maps, by considering the average distance of all paths identified by the software. The distance matrices are presented in the Appendix, where the suppliers, customers, and distribution centres are indicated by letters f, c, and d, respectively, for simplicity. In detail, f1 = Marsala, f2 = Menfi, f3 = Bronte, f4 = Modica, f5 = Capo d’Orlando, f6 = Piazza Armerina, f7 = Riposto, d1 = Alcamo, d2 = Acireale, d3 = Vittoria, d4 = San Cataldo, c1 = Palermo, c2 = Messina, c3 = Trapani, c4 = Catania, c5 = Agrigento, c6 = Enna, c7 = Caltanissetta, c8 = Siracusa, and c9 = Ragusa.

Journal Pre-proof

Figure 2: Geographical representation of supply chain actors Table 1: Sets and indices of 3E-FGVRP i

Index of customer groups: i ϵ I = {1, 2,…, c,…, C}, where C is the number of customers.

j

Index of food suppliers: j ϵ J = {1, 2,…, f,…, F}, where F is the number of suppliers.

k

Index of distribution centres: k ϵ K = {1, 2,…, d,…, D}, where D is the number of depots.

n

Index of vehicles: n ϵ N = {1, 2,…, v,…, V}, where V is the vehicle fleet size.

t

Index of time periods: t ϵ T = {1, 2,…, t,…, T} 𝑇,where T is the number of time periods.

A

E

U

H

L

Set of arcs connecting every pair of suppliers j, j ϵ J. A = {a1, a2,…, a−1,…, a}𝐴 =

{𝑎1,𝑎2,…,𝑎𝐴 ― 1,𝑎𝐴}, where 𝐴 is the number of arcs (j, j). Set of arcs connecting suppliers and distribution centres j∈ J, k∈ K. 𝐸 = {𝑒1,𝑒2,…,𝑒𝐸 ― 1,𝑒𝐸}, where 𝐸 is the number of arcs (j, k). Set of arcs connecting every pair of distribution centres k, k ∈ K. 𝑈 = {𝑢1,𝑢2,…,𝑢𝑈 ― 1,𝑢𝑈}, where 𝑈 is the number of arcs (k, k). Set of arcs connecting every pair of customer groups i, i ∈ I. 𝐻 = {ℎ1,ℎ2,…,ℎ𝐻 ― 1,ℎ𝐻}, where 𝐻 is the number of arcs (i, i). Set of arcs connecting distribution centres and customer groups k∈ K, i∈ I. 𝐿 =

{𝑙1,𝑙2,…,𝑙𝐿 ― 1,𝑙𝐿}, where 𝐿 is the number of arcs (k, i).

Z Set of all supply chain arcs: 𝑍 = 𝐴 ∪ 𝐸 ∪ 𝑈 ∪ 𝐻 ∪ 𝐿 indexed by z.

Journal Pre-proof Table 2: Parameters Symbol

Explanation

Unit

Mk

Capacity of k-th distribution centre.

[items]

Qn

Capacity of n-th transportation vehicle.

[items]

Total number of food products requested by i-th customer from all

[items]

𝑑𝑒𝑚𝑡𝑖 𝐷𝑒𝑚𝑡𝑗

suppliers at period t. Set of food products requested by all customer groups from j-th supplier.

[items]

FCk

Fixed cost of opening a distribution centre k.

[€]

Fcn

Fixed cost of employing a vehicle n.

[€]

TCz

Travelling costs associated with arc z ∈ Z.

[€/km]

δz

Travel distance associated with arc z ∈ Z.

[km]

Ef

CO2 emissions when vehicle is full.

[kg/km]

Ee

CO2 emissions when vehicle is empty.

[kg/km]

𝛼,𝛽,𝛾

Preference thresholds of fuzzy chance constraints.

Before presenting the fuzzy chance-constrained programming model, the main assumptions related to the 3E-FGVRP mathematical model are given below. 

Both vehicle and distribution centres can store a limited amount of goods; therefore, the total load must not exceed their physical capacity.



A homogeneous and limited fleet of vehicles is available at each period to pick up and deliver food products.



A vehicle is assigned for only one route at period t; therefore, at the end of each period, the entire fleet is completely available to deliver other foodstuffs at t+1.



The suppliers and customers are visited by one and only one vehicle at each time period.



All vehicles are initially empty and start the routes from a distribution centre by picking up the total amount of goods requested by all customers from a supplier. After visiting a certain number of suppliers (for a single route), the vehicles return to a depot.



The overall demand of each customer must be assigned to the same distribution centre; this is the reason that additional sub-routes between depots are planned.



In the third echelon of the supply chain, every route starts in distribution centres, continues by visiting a certain number of customers, and finishes in a distribution centre.



The itineraries among food suppliers, distribution centres, and customer groups are represented by arches to which the real distance measures are associated.

Journal Pre-proof 

The demand of a customer at each period is a triangular fuzzy number, represented by the following notation: 𝑑𝑒𝑚𝑡𝑖 = (𝑑𝑒𝑚1,𝑖𝑡,𝑑𝑒𝑚2,𝑖𝑡,𝑑𝑒𝑚3.𝑖𝑡). The detailed description of the customer demand simulation algorithm is presented in Section 4.

Although the capacity of each vehicle Qn is a deterministic quantity, it is not possible to plan the routes a priori. This happens because in the design preliminary phase, only an approximate information of the product amount (not the exact quantity) requested by the customers is known. After having picked up the overall demands from supplier f, the residual capacity of the n-th vehicle 𝑓

is 𝑅𝑡𝑛,𝑓 = 𝑄𝑛 ― ∑𝑗 = 1𝐷𝑒𝑚𝑡𝑗 for the routes of the supply chain’s first echelon. Obviously, 𝐷𝑒𝑚𝑡𝑗 is also a fuzzy triangular number as this is the sum of all C customer demands related to the food product provided by supplier f. These demands and the residual capacity of each vehicle are obtained by using the basic operations on fuzzy relations [84]. In particular, Eq. (5) shows the extended form of the fuzzy residual capacity of the n-th vehicle after visiting f suppliers; for simplicity in the discussion, Eq. (5) is represented by a shorter form, i.e. Eq. (6).

(

𝑓

𝑓

𝑓

𝑅𝑡𝑛,𝑓 = 𝑄𝑛 ― ∑𝑗 = 1𝐷𝑒𝑚𝑡3,𝑗 ,𝑄𝑛 ― ∑𝑗 = 1𝐷𝑒𝑚𝑡2,𝑗,𝑄𝑛 ― ∑𝑗 = 1𝐷𝑒𝑚𝑡1,𝑗

)

𝑅𝑡𝑛,𝑓 = (𝑟1,𝑛𝑓𝑡,𝑟2,𝑓𝑡,𝑟3.𝑓𝑡)

(5) (6)

In light of the above, the fuzzy event ‘pick up from next supplier’ is associated with a probability of occurrence that is calculated by means of the credibility theory as already mentioned. The logical steps necessary for computing the credibility measure related to the first echelon of supply chain are given by Eq. (7). A discussion concerning the first stage of the supply chain is given below, and the same expressions and considerations are valid for the third echelon, with the difference that the residual capacity of each vehicle is calculated by considering the overall foodstuffs requested by each 𝑐

customer 𝑅𝑡𝑛,𝑐 = 𝑄 ― ∑𝑖 = 1𝑑𝑒𝑚𝑡𝑖 . A credibility measure for the fuzzy event ‘serve next customer’ has to be found.

{

}

𝐶𝑟𝑒𝑑 = 𝐶𝑟 𝐷𝑒𝑚𝑡𝑓 + 1 ≤ 𝑅𝑡𝑛,𝑓

= 𝐶𝑟{(𝐷𝑒𝑚𝑡1,𝑓 + 1 ― 𝑟3,𝑓𝑡,𝐷𝑒𝑚𝑡2,𝑓 + 1 ― 𝑟2,𝑓𝑡,𝐷𝑒𝑚𝑡3,𝑓 + 1 ― 𝑟1,𝑓𝑡) ≤ 0}

{

𝑟3,𝑓𝑡

0 ― 𝐷𝑒𝑚𝑡1,𝑓 + 1

2(𝑟3,𝑓𝑡 ― 𝐷𝑒𝑚𝑡1,𝑓 + 1 + 𝐷𝑒𝑚𝑡2,𝑓 + 1 ― 𝑟2,𝑓𝑡) = 𝐷𝑒𝑚𝑡3,𝑓 + 1 ― 𝑟1,𝑓𝑡 ― 2(𝐷𝑒𝑚𝑡2,𝑓 + 1 ― 𝑟2,𝑓𝑡)

𝑖𝑓 𝐷𝑒𝑚𝑡1,𝑓 + 1 ≥ 𝑟3,𝑓𝑡 𝑖𝑓 𝐷𝑒𝑚𝑡1,𝑓 + 1 ≤ 𝑟3,𝑓𝑡,𝐷𝑒𝑚𝑡2,𝑓 + 1 ≥ 𝑟2,𝑓𝑡

𝑖𝑓 𝐷𝑒𝑚𝑡2,𝑓 + 1 ≤ 𝑟2,𝑓𝑡,𝐷𝑒𝑚𝑡3,𝑓 + 1 ≥ 𝑟1,𝑓𝑡 2(𝑟2,𝑓𝑡 ―𝐷𝑒𝑚𝑡2,𝑓 + 1 + 𝐷𝑒𝑚𝑡3,𝑓 + 1 ― 𝑟1,𝑓𝑡) 1 𝑖𝑓 𝐷𝑒𝑚𝑡3,𝑓 + 1 ≤ 𝑟1,𝑓𝑡

(7)

Journal Pre-proof

Similarly, for the vehicle capacity, having considered the customer demands as uncertain data, the residual capacity of the k-th distribution centre at period t after allocating the overall demands of c customer is also a triangular fuzzy number presented in Eq. (8).

(

𝑐

𝑐

𝑐

)

𝑃𝑡𝑘,𝑐 = 𝑀𝑘 ― ∑𝑖 = 1𝑑𝑒𝑚𝑡3,𝑖 ,𝑀𝑘 ― ∑𝑖 = 1𝑑𝑒𝑚𝑡2,𝑖,𝑀𝑘 ― ∑𝑖 = 1𝑑𝑒𝑚𝑡1,𝑖

(8)

The latter expression is abbreviated below with a simpler form, i.e. Eq. (9), in which the credibility measure corresponding to the event ‘next customer’s demand is less than or equal to the depot’s available capacity’ is calculated by using Eq. (10). 𝑃𝑡𝑘,𝑐 = (𝑝1,𝑘𝑡,𝑝2,𝑘𝑡,𝑝3.𝑘𝑡)

{

(9)

}

𝐶𝑟𝑒𝑑 = 𝐶𝑟 𝑑𝑒𝑚𝑡𝑐 + 1 ≤ 𝑃𝑡𝑘,𝑐

(10)

= 𝐶𝑟{(𝑑𝑒𝑚𝑡1,𝑐 + 1 ― 𝑝3,𝑘𝑡,𝑑𝑒𝑚𝑡2,𝑐 + 1 ― 𝑝2,𝑘𝑡,𝑑𝑒𝑚𝑡3,𝑐 + 1 ― 𝑝1,𝑘𝑡) ≤ 0}

{

𝑝3,𝑘𝑡

0 ― 𝑑𝑒𝑚𝑡1,𝑐 + 1

2(𝑝3,𝑘𝑡 ― 𝑑𝑒𝑚𝑡1,𝑐 + 1 + 𝑑𝑒𝑚𝑡2,𝑐 + 1 ― 𝑝2,𝑘𝑡) = 𝑑𝑒𝑚𝑡3,𝑐 + 1 ― 𝑝1,𝑘𝑡 ― 2(𝑑𝑒𝑚𝑡2,𝑐 + 1 ― 𝑝2,𝑘𝑡)

𝑖𝑓 𝑑𝑒𝑚𝑡1,𝑐 + 1 ≥ 𝑝3,𝑘𝑡 𝑖𝑓 𝑑𝑒𝑚𝑡1,𝑐 + 1 ≤ 𝑝3,𝑘𝑡,𝑑𝑒𝑚𝑡2,𝑐 + 1 ≥ 𝑝2,𝑘𝑡

𝑖𝑓 𝑑𝑒𝑚𝑡2,𝑐 + 1 ≤ 𝑝2,𝑘𝑡,𝑑𝑒𝑚𝑡3,𝑐 + 1 ≥ 𝑝1,𝑘𝑡 2(𝑑𝑒𝑚𝑡2,𝑐 + 1 ― 𝑝2,𝑘𝑡 + 𝑑𝑒𝑚𝑡3,𝑐 + 1 ― 𝑝1,𝑘𝑡) 1 𝑖𝑓 𝑑𝑒𝑚𝑡3,𝑓 + 1 ≤ 𝑝1,𝑘𝑡

The process of determining routes is based on very intuitive considerations: if the quantity picked up from the j-th supplier is fairly moderate and the overall demand of one of the suppliers not yet visited is small enough, the credibility that the vehicle has sufficient space to deposit this quantity is relatively high (very close to 1). On the contrary, the smaller the difference between the vehicle’s residual available capacity and demand of the next suppliers, the smaller the credibility value (most likely close to 0), i.e. the probability of sending the vehicle to the next supplier is very small. The cases above represent the limit conditions under which the vehicle will certainly be sent to the next supplier or will return to the depot. Therefore, some intermediate cases characterised by Cred ∈ (0, 1) can occur. In these circumstances, the question of whether the vehicle will continue the route or return to a distribution centre arises. Thus, it is necessary to introduce a preference threshold α representing the decisive limit credibility for route design. The parameter α is the upper bound of a confidence interval containing the plausible values for the credibility measure. It should be noted that the preference threshold α is chosen subjectively by the supply chain manager in the [0, 1] range as it represents the manager’s attitude towards risk. The closer the value to 1, the lower the risk of going

Journal Pre-proof to the next supplier with an insufficient capacity. In this case, the routes are defined by adopting a prudential policy, characterised by a less utilisation of vehicle capacity. On the other hand, if the supply chain manager chooses an α value close to 0, the manager pursues a risky policy, meaning that he/she wants to maximise the quantity picked up at each route and minimise the total number of routes. In the latter case, the supply chain manager plans routes by using the vehicle’s available capacity as much as possible. Thus, according to the above considerations, the condition Cred ≥ α is critical for deciding whether to send the vehicle to the next supplier or to return it to a distribution centre. If Cred ≥ α holds, the vehicle is sent to the next supplier; otherwise, it is brought back to a distribution centre and a different fleet of vehicle is sent to the supplier. The mathematical model is set up by repeating the route definition until all the suppliers are visited by a vehicle. Moreover, for the third echelon of the supply chain, the preference threshold γ is defined. The same consideration can be made on the credibility measures for distribution centres. When Cred = 0, the supply chain manager is aware that the distribution centre k at period t is certainly not able to satisfy the demand of another customer. In contrast, when Cred = 1, the depot is able to satisfy the quantity. Once again, it is necessary to set up a preference threshold to determine whether to stock the next customer’s demand in depot k or not. If Cred ≥ β is fulfilled, it is credible that the depot’s available capacity is less than or equal to the next customer’s demand, and therefore, the distribution centre k at period t should store the quantity; otherwise, the customer should be served by another open depot. The values of preference thresholds involved in this study are provided in Section 7 and are set up from the contributions of sector experts. The selected threshold values represent any trade-off solutions between costs and carbon emissions and customer satisfaction. Being a B2C supply chain, it is crucial to plan the service in the most efficient way to ensure that customers are highly satisfied by suggesting other customers to use the service. For the aforementioned reason, a less risky policy is adopted. 3.1 Mathematical formulation of 3E-FGVRP This section presents the mathematical model of the 3E-FGVRP by describing the decision variables, objective functions, and constraint equations of the corresponding chance-constrained programming model based on the credibility theory. 3.1.1 Decision variables In accordance with the assumptions above, the following binary decision variables are defined.

𝑋𝑡𝑧,𝑛 =

𝑡 𝑡𝑟𝑎𝑣𝑒𝑙𝑠 𝑜𝑛 𝑎𝑟𝑐ℎ 𝑧 ∈ {𝑍} {1 𝑖𝑓 𝑣𝑒ℎ𝑖𝑐𝑙𝑒 𝑛 𝑎𝑡 𝑝𝑒𝑟𝑖𝑜𝑑 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

Journal Pre-proof 𝑘 𝑎𝑡 𝑝𝑒𝑟𝑖𝑜𝑑 𝑡 𝑖𝑠 𝑜𝑝𝑒𝑛 {1 𝑖𝑓 𝑑𝑖𝑠𝑡𝑟𝑖𝑏𝑢𝑡𝑖𝑜𝑛0,𝑐𝑒𝑛𝑡𝑟𝑒 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 1 𝑖𝑓 𝑑𝑒𝑚𝑎𝑛𝑑 𝑜𝑓 𝑐𝑢𝑠𝑡𝑜𝑚𝑒𝑟 𝑖 𝑎𝑡 𝑝𝑒𝑟𝑖𝑜𝑑 𝑡 𝑖𝑠 𝑠𝑒𝑟𝑣𝑒𝑑 𝑏𝑦 𝑑𝑖𝑠𝑡𝑟𝑖𝑏𝑢𝑡𝑖𝑜𝑛 𝑐𝑒𝑛𝑡𝑟𝑒 𝑘 Γ ={ 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 1 𝑖𝑓 𝑣𝑒ℎ𝑖𝑐𝑙𝑒 𝑛 𝑎𝑡 𝑝𝑒𝑟𝑖𝑜𝑑 𝑡 𝑖𝑠 𝑠𝑒𝑙𝑒𝑐𝑡𝑒𝑑 ζ ={ 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 Υ𝑡𝑘 = 𝑡 𝑖,𝑘

𝑡 𝑛

3.1.2 Objective functions In line with the classic GVRP applications, this study aims at designing a distribution supply chain according to two objective functions, namely total costs and CO2 emissions, which must be minimised. Both formulations are provided below in Eqs. (11) and (12), respectively. 3.1.2.1 Cost function In this study, the overall cost function consists of the sum of two fixed cost items and the last is related to variable costs. In particular, the fixed costs regard the expenditure involved in employing distribution centres and vehicles. The parameters FCk and Fcn are calculated by considering all possible variable and fixed costs for the k-th distribution centre and n-th vehicle (e.g. purchase costs of vehicles, cost of utilities, employment of personnel, building rental for the distribution centre, and costs of maintenance and personnel for fleet vehicles). The transportation costs are variable because they are directly proportional to the distance characterising every travelled arch. The parameter TC represents the operating costs for a vehicle that approximately covers 100,000 km per year. This parameter, which includes the average cost of fuel, corresponds to the real data of the Italian automobile club. The mathematical formulation of the total cost function is given by Eq. (11) [85]. 𝑀𝑖𝑛𝑖𝑚𝑖𝑧𝑒: 𝑇𝑜𝑡𝑎𝑙𝐶𝑜𝑠𝑡 = ∑𝑡 ∈ 𝑇∑𝑘 ∈ 𝐾𝐹𝐶𝑘Υ𝑡𝑘 + ∑𝑡 ∈ 𝑇∑𝑛 ∈ 𝑁𝐹𝑐𝑛ζ𝑡𝑛 + ∑𝑡 ∈ 𝑇∑𝑛 ∈ 𝑁∑𝑧 ∈ 𝑍𝑇𝐶 𝛿𝑧𝑋𝑡𝑧,𝑛

(11)

3.1.2.2 Carbon emission function The formulation of the carbon emission function in this study is derived from a careful analysis of the literature. In particular, the amount of carbon dioxide released by the vehicle fleet is calculated by considering references [86–88], although with some modifications. Those studies considered the freight transport for heavy-duty vehicles with a capacity between 32 and 40 tons under some important assumptions, such as a constant average vehicle speed of 80 km/h and a perfectly flat road. In this study, the overall expected transport volumes are clearly less than the ones in the aforementioned references; thus, the parameters Ef and Ee (which represent the carbon emissions under the conditions of full and empty load, respectively) are computed by customising them to the 3E-FGVRP. In addition, the general formula for calculating the emission at each arc and the

Journal Pre-proof assumptions are considered exactly in the same way. Eq. (12) provides the mathematical expression by which the CO2 emissions are evaluated in terms of the distance (km) and weight of quantity loaded (kg) at each arc. Because this study considers the quantities requested by customers as fuzzy data, 𝑞𝑧𝑡 and 𝑄𝑧,𝑛𝑡 are considered as the medium points of triangular fuzzy numbers corresponding to food products loaded in a generic arc z at period t and the capacity of vehicle n at period t and arc z, respectively.

(

𝑞𝑧𝑡

)

Minimise: CarbonEmission = ∑𝑡 ∈ 𝑇∑𝑛 ∈ 𝑁∑𝑧 ∈ 𝑍𝛿𝑧𝑋𝑡𝑧,𝑛 (𝐸𝑓 ― 𝐸𝑒)𝑄 𝑡 + 𝐸𝑒 𝑧,𝑛

(12)

3.1.3 Problem constraints The constraint equations to which the objective functions are subjected are presented below. ∑𝑛 ∈ 𝑁𝑋𝑡𝑧,𝑛 = 1 ∀ 𝑡 ∈ 𝑇, ∀ 𝑧 ∈ 𝐴, 𝐸, 𝐻, 𝐿 ⊆ 𝑍

(13)

∑𝑧 ∈ 𝐸 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ― ∑𝑧 ∈ 𝐸 ⊆ 𝑍𝑋𝑡𝑧,𝑛 = 0

∀ 𝑡 ∈ 𝑇 ,∀ 𝑛 ∈ 𝑁

(14)

∑𝑧 ∈ 𝐿 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ― ∑𝑧 ∈ 𝐿 ⊆ 𝑍𝑋𝑡𝑧,𝑛 = 0

∀ 𝑡 ∈ 𝑇 ,∀ 𝑛 ∈ 𝑁

(15)

( 𝐶𝑟𝑒𝑑(∑ 𝐶𝑟𝑒𝑑(∑

) ≤𝑄 )≥𝛾

𝐶𝑟𝑒𝑑 ∑𝑗 ∈ 𝐽∑𝑧 ∈ 𝐴 ⊆ 𝑍𝐷𝑒𝑚𝑡𝑗𝑋𝑡𝑧,𝑛 ≤ 𝑄𝑛 ≥ 𝛼

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(16)

∑𝑧 ∈ 𝐻 ⊆ 𝑍𝑑𝑒𝑚𝑡𝑖𝑋𝑡𝑧,𝑛

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(17)

𝑖∈𝐼

𝑛

)

𝑑𝑒𝑚𝑡𝑖Γ𝑡𝑖,𝑘 ≤ 𝑀𝑘Υ𝑡𝑘 ≥ 𝛽

𝑖∈𝐼

∀ 𝑡 ∈ 𝑇, ∀ 𝑘 ∈ 𝐾

(18)

∑𝑛 ∈ 𝑁∑𝑧 ∈ 𝐴 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝐴

∀𝑡∈𝑇

(19)

∑𝑛 ∈ 𝑁∑𝑧 ∈ 𝐻 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝐻

∀𝑡∈𝑇

(20)

∑𝑛 ∈ 𝑁𝑋𝑡𝑧,𝑛 ≤ 𝑉 ∀ 𝑡 ∈ 𝑇, ∀ 𝑧 ∈ 𝑈 ⊆ 𝑍

(21)

∑𝑧 ∈ 𝑈 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝑈 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(22)

∑𝑧 ∈ 𝐴 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝐹 ― 1 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(23)

∑𝑧 ∈ 𝐻 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝐶 ― 1 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(24)

∑𝑧 ∈ 𝐴

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

𝑡

1

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁,∀ 𝐴 ∈ 𝐴 ⊆ 𝑍

(25)

∑𝑧 ∈ 𝐻

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

𝑡

1

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁,∀ 𝐻 ∈ 𝐻 ⊆ 𝑍

(26)

∑𝑧 ∈ 𝐸

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

𝑡

1 𝑎𝑛𝑑 ∑𝑧 ∈ 𝐸

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

∑𝑧 ∈ 𝐿

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

𝑡

1 𝑎𝑛𝑑 ∑𝑧 ∈ 𝐿

𝑋 ≤ ⊆ 𝑍 𝑧,𝑛

∑𝑖 ∈ 𝐼Γ𝑡𝑖,𝑘 ― ∑𝑧 ∈ 𝑈 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 𝐶 ― 1

𝑡

1 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁,∀ 𝐸 ∈ 𝐸 ⊆ 𝑍, ∀ 𝐸 ∈ 𝐸 ⊆ 𝑍

(27)

𝑡

1 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁, ∀ 𝐿 ∈ 𝐿 ⊆ 𝑍 , ∀ 𝐿 ∈ 𝐿 ⊆ 𝑍

(28)

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁, ∀ 𝑘 ∈ 𝐾

(29)

Journal Pre-proof ∑𝑧 ∈ 𝐸 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ― ∑𝑧 ∈ 𝑈 ⊆ 𝑍𝑋𝑡𝑧,𝑛 ≤ 1 𝑋𝑡𝑧,𝑛 ∈ {0,1} Υ𝑡𝑘 ∈ {0,1} Γ𝑡𝑖,𝑘 ∈ {0,1} ζ𝑡𝑛 ∈ {0,1}

(30)

∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁, ∀ 𝑘 ∈ 𝐾

(31)

∀ 𝑡 ∈ 𝑇, ∀ 𝑧 ∈ 𝑍, ∀ 𝑛 ∈ 𝑁 ∀ 𝑡 ∈ 𝑇, ∀ 𝑘 ∈ 𝐾

(32) (33)

∀ 𝑡 ∈ 𝑇, ∀ 𝑘 ∈ 𝐾, ∀ 𝑖 ∈ 𝐼 ∀ 𝑡 ∈ 𝑇, ∀ 𝑛 ∈ 𝑁

(34)

Constraint (13) ensures that an arc 𝑧 ∈ 𝐴, 𝐸, 𝐻, 𝐿 ⊆ 𝑍 must be visited by only one vehicle during the same period t. Eqs. (14) and (15) guarantee the continuity of routes for the vehicles that travel on paths belonging to subsets 𝐸,𝐿 , 𝐸, and 𝐿, representing the arcs travelled during the outward and the return journey, respectively. The fuzzy chance constraints are presented in Eqs. (16)–(18). Eqs. (16)– (17) ensure that all suppliers and customers are visited within vehicle capacity and Eq. (18) guarantees that each customer’s overall demand is allocated within the distribution centre capacity with defined confidence levels. Constraints (19) and (20) ensure that the generic arc z is visited once in every period t at the most. Eq. (21) ensures that an arc 𝑧 ∈ 𝑈 ⊆ 𝑍 can be travelled by more than one vehicle (at most V) in the same period t, because U represents the set of additional arcs between distribution centres. Constraints (22)–(24) define the maximum number of arcs that a vehicle can visit in every period t for different arc sets, respectively. Eqs. (25)–(28) ensure that a vehicle travels only one arc among those related to the same supply chain actor in period t (these subsets are represented by the double line accent, e.g. 𝐴). The relations indicated in Eqs. (29)–(30) ensure the continuity of routes of each vehicle n at period t. Finally, constraints (31)–(34) specify the binary variables used in the formulation. 4. Customer demand simulation This study aims at solving a 3E-FGVRP by applying a heuristic multi-objective optimisation algorithm in determining the efficient routes of vehicles for a three-stage supply chain in the Sicilian agri-food sector. The input data necessary to conduct the study are the quantities of products requested by each customer group. As a matter of fact, in any real context, an order receipt represents the starting point from which the distribution dynamics of the supply chain begin. As this is a preliminary study, the data related to the amount of foodstuffs requested by customers are simulated by taking into consideration some details that represent the real consumers’ behaviour for purchasing food. In particular, the authors decided to include the uncertainty associated with the variability of the required quantities by generating customer demands as fuzzy triangular numbers. To the best of the authors’ knowledge, this is the first study in which a GVRP is solved by using nondeterministic data derived from a simulation algorithm.

Journal Pre-proof First, the code is applied to simulate the real quantities of products requested by customers at each period t. Then, the evolutionary algorithm is executed to solve the above fuzzy chance-constrained problem to which vehicles’ routes according to the supply chain minimum costs and CO2 emissions are found. As mentioned above, this study considers several consecutive periods during which the route definition is carried out independently of each other (i.e. the receipt of orders and related distributions must be carried out entirely during a single period). However, these simulations are highly dependent on one another due to consumers’ purchasing behaviour, that is, the amount of goods requested by the consumer in different periods. The implemented simulation algorithm takes into consideration the number of products that the customer orders from supplier f at time t−1. To re-create the consumer’s purchasing behaviour as closely as possible to reality, the quantities of product ordered by a consumer at time t are proportional to the same goods quantity ordered in the previous period. Therefore, if customer c has ordered a relatively high number of products at time t−1 from supplier f, most likely, at period t, the same customer will order a smaller number of products. Similarly, if at time t−1 the customer orders the quantity y of products, he/she requires the quantity g at period t (where g < y) and at t+1, he/she orders the quantity w (w > g). The above assertions refer to the median values of the triangular fuzzy numbers, i.e. those of the triangle with the highest probability of occurrence. Therefore, the following study allows designing in advance the supply chain (number of vehicles employed and relative capacity; size of distribution centres). In fact, it is well known that the values of demand and relative temporal variation become certain a posteriori when the customer places the order, being completely sure of the quantities he/she wants to order and consume. To take into consideration the uncertainty and variability that commonly characterise the purchasing behaviour of consumers, the implemented simulation algorithm returns fuzzy demand values, as shown below. The customers demand simulation is summarised as follows. Step 1: The behavioural profiles of the customers (in terms of the number of products purchased from each supplier) are identified. In particular, these key data are extrapolated by an introductory analysis of the agri-food sector concerning the benchmarking of the main Italian e-commerce platforms and cooperation with several sector experts. These collaborations and studies lead to the definition of the demand simulation’s input data, i.e. the ranges containing all possible values of requested foodstuffs by a customer for each supplier. From the latter, the values representing each fuzzy triangle number (lower, medium, and upper) are randomly extrapolated.

Journal Pre-proof Step 2: The procedure of Step 1 returns the demand’s fuzzy triangular numbers for the first period t. By following the guidelines mentioned above, the code implemented in Step 2 defines iteratively the demand values related to the next period by considering the total amount of previous orders. Hence, the temporal deviations of the demands of the subsequent periods are calculated from iteration to iteration. 5. NSGA-II for 3E-FGVRP resolution NSGA-II is one of the most popular multi-objective evolutionary algorithms (MOEAs) and it has been applied in many engineering applications. The main feature of NSGA-II is the elitist sorting method in which the crowding distance is used as the ranking criterion for sorting the non-dominated solutions. 5.1 Description NSGA-II is an optimisation algorithm that belongs to the evolutionary algorithm family. It proposes an extension of the genetic algorithm for multiple objective function optimization. The NSGA aims at improving the adaptive fit of a population of candidate solutions to a Pareto front constrained by a set of objective functions [89]. The algorithm uses an evolutionary process made up of different evolutionary operators including selection, genetic crossover, and genetic mutation. In the first step, the population is sorted into several sub-populations based on Pareto dominance. Successively, the similarity between members of each Pareto sub-group is evaluated by means of a specific measure, which is used to promote a diverse front of non-dominated solutions. In the second step, the selected parents are combined using the crossover operator to create the offspring chromosomes. The latter are finally mutated to guarantee that the population is composed of individuals with a certain degree of diversity. 5.1.1 Coding of initial population The genetic algorithm is a search method that requires as input data a population of feasible solutions, named the initial population, rather than a single one. A specific number of solutions are randomly created, being this latter the pre-defined population size Pop_size. Moreover, all solutions belonging to the initial population have to satisfy the declared problem constraints. Each generation of genetic minimisation produces a population of Pop_size feasible solutions; therefore, this size does not change during the optimisation. However, the procedure described below is necessary only for the creation of initial population.

Journal Pre-proof Each individual in the population (also named chromosome in the technical language of genetic algorithms) consists of θ genes representing the number of problem decision variables (the value θ will be analogously defined as genome length). In particular, a chromosome is represented by a row vector of dimension (1, θ), whose elements are uniquely assigned to the value of a specific decision variable. As previously mentioned, in this study, all decision variables are set as binary; thus, the implemented random code by which the initial population chromosomes are developed is as follows. Chromosome (row,:) = round (rand (1, θ)), row = 1, 2,…, Pop_size

(35)

Before it can be accepted as a possible solution, it is necessary to verify that the generated chromosome meets the constraints of the problem. If the constraints are not violated, the chromosome is accepted and stored as a feasible solution of the initial population array; otherwise, the creation code related to the unfeasible chromosome is repeated until all the constraints are verified. 5.1.2 Elitist ranking chromosomes Once the initial population is created, the individuals are sorted by following the non-domination principle. The chromosome ranking is performed using elitism, which helps to achieve better convergence in MOEAs [90]. Moreover, elitism allows to carry over the best chromosomes unchanged, from the current population to the next ones, by guaranteeing that the quality of solutions will not decrease from one generation to the next. The elitist procedure is combined with selection to ensure that in every generation the best chromosomes are more likely to be selected. Elitism consists of two phases: the aggregation of chromosomes into classes according to the nondomination principle, and the chromosomes ranking in agreement with both the importance of class to which they belong to and their measure of proximity (similarity) with the others. 5.1.2.1 Non-dominated sorting An individual dominates another one if all of its objective function values are no worse than the other ones and if at least one of its objective functions is better than the other ones. According to the statement above, corresponding to the non-dominance principle [91], each chromosome is compared with all the other ones to verify whether this is a dominant solution or not. All dominating individuals not dominated by anyone else are assigned to a class corresponding to Pareto’s first frontier. This classification procedure is repeated iteratively until all the solutions have been classified. Step by step, the number of chromosomes to be compared becomes increasingly smaller because the solutions already assigned to a Pareto front are not considered in the next iterations of classification. Moreover,

Journal Pre-proof a rank value is associated with each frontier, being Rank = 1 for the first front, Rank = 2 for the second, and Rank = x for the x-th defined Pareto front. This procedure allows creating a first rough ranking of chromosomes because the solutions characterised by a lower rank value are considered better than the others. Therefore, all the individuals of the first frontier are better than those of the second, which in turn are better than the individuals of the third front and so on. At this phase, nothing can be said about the solutions’ quality belonging to the same frontier. 5.1.2.2 Crowding distance The overall ranking of the Pop_size chromosomes of the population is defined through the crowding distance that allows sorting the solutions characterised by the same rank. The crowding distance evaluates the density degree of a solution compared to that of other solutions belonging to the same front. In particular, the distance computation is performed by considering each objective function at a time and eventually aggregating the results. Therefore, first, the chromosomes of the same front are sorted in ascending order according to the value of each objective function. Subsequently, an infinite distance value is assigned to the boundary individuals (i.e. chromosomes with the smallest and the largest function values), whereas the Euclidean distance is used to compute the distance between two adjacent individuals. In particular, because the optimisation problem analysed in this study has two objective functions, the crowding distance is the diagonal length of a rectangle formed by using the adjacent individuals as the vertices. The analytical approach of the steps just described are reported in [92]. The i-th chromosome of population is sorted by means of the crowded operator ( ≺ ), which takes into consideration both non-dominated rank (irank) and crowding distance (icrowding). Eq. (36) summarises the rules by which the ranking of individuals is determined. Chromosome i is preferred to chromosome j if its non-dominated rank is lower than that of j (i.e. i belongs to a better Pareto front) or if both the solutions belong to the same Pareto front and the crowding distance of chromosome i is greater than that of j (i.e. a solution located in a less crowded region is preferred). 𝑖 ≺ 𝑗 𝑖𝑓 (𝑖𝑟𝑎𝑛𝑘 < 𝑗𝑟𝑎𝑛𝑘) ∨ (𝑖𝑟𝑎𝑛𝑘 = 𝑗𝑟𝑎𝑛𝑘 ∧ 𝑖𝑐𝑟𝑜𝑤𝑑𝑖𝑛𝑔 > 𝑗𝑐𝑟𝑜𝑤𝑑𝑖𝑛𝑔 )

(36)

5.2 Selection Once the ranking is defined, the parents must be chosen to create the offspring population. The chromosome selection method applied in this study is the roulette wheel with some elements of novelty. The procedure for the selection of an individual’s probability measure is described below.

Journal Pre-proof First, a descending parameter named rank_front is associated with each Pareto frontier as follows: if x front is found, the frontier characterised by Rank = 1 has a rank_front = x, a rank_front = x−1 is associated with the Pareto front with Rank = 2, etc. All chromosomes belonging to the same front are characterised by a single rank_front value. Second, the procedure just described is repeated internally to individuals on each front. If the generic x-th front is composed of w chromosomes, the parameter rank_chromosome = w is assigned to the individual ranked first, a rank_chromosome = w−1 is associated with the second best chromosome, and so on, until the rank_chromosome = 1 is given to the worst chromosome of the front. By summing up the rank_front and rank_chromosome of every solution and by normalising them with respect to the total, the selection probabilities are obtained. By doing so, the best solutions are more likely to be selected, but they will not necessarily be chosen. As a matter of fact, similar to a roulette game, the parent selection is associated with a random choice to avoid the situation in which during the next iterations, the same individuals might not often be chosen. After generating a random number in [0; 1] from time to time, the selected parent is the one whose probability of selection is the closest to its value. Because this random procedure could repeatedly select the same parent chromosomes, a check step is inserted in the algorithm code: if a parent has already been selected in the same iteration, it will no longer be chosen. Therefore, another random number is generated to select a different parent. Figures 3–4 and Table 3 present a representative example of how parent chromosomes are selected. In particular, a sample population of 21 feasible individuals is considered. These individuals are previously classified into four Pareto fronts and successively ranked by means of non-dominated sorting and crowding distance. Figure 3 shows the scatter plot of the individuals with regard to the objective function values. As already mentioned, the chromosomes belonging to the same frontier (represented with the same markers and colours) are characterised by the rank_front value indicated in the figure legend and by a rank_chromosome, whose values are shown above each symbol.

Journal Pre-proof

Figure 3: Pareto fronts of sample population chromosomes Finally, Table 3 presents the method in which the selection probability associated with each solution is calculated, whereas Figure 4 shows the roulette, whose section corresponds to the selection probability of each chromosome. Table 3: Computation of individuals’ selection probabilities Rank front Rank

Rank

Normalised rank

front

chromosome

chromosome

+ Normalised rank

Selection probability

chromosome 1st Front

2nd Front

3rd Front

4

3

2

4

1

5

0.0813

3

0.75

4.75

0.0772

2

0.5

4.5

0.0732

1

0.25

4.25

0.0691

5

1

4

0.0650

4

0.8

3.8

0.0618

3

0.6

3.6

0.0585

2

0.4

3.4

0.0553

1

0.2

3.2

0.0520

6

1

3

0.0488

5

0.83

2.83

0.0461

4

0.67

2.67

0.0434

3

0.5

2.5

0.0407

2

0.33

2.33

0.0379

1

0.17

2.17

0.0352

Journal Pre-proof

4th Front

1

6

1

2

0.0325

5

0.83

1.83

0.0298

4

0.67

1.67

0.0271

3

0.5

1.5

0.0244

2

0.33

1.33

0.0217

1

0.17

1.17

0.0190

61.5

Figure 4: Representation of probability roulette 5.3 Crossover The genetic evolutionary search is mainly based on a crossover operator that allows generating the offspring population. In the literature, several ways of performing the crossover operator exist [93]. In this study, the authors decide to implement an x-point crossover, where x is a fraction of the genome length. The points where parent strings are involved during the crossover are determined randomly by creating a row vector (1, x), named crossover points of integer random numbers in the range [1, genome length]. Figure 5 shows the method of creating the two new strings. By taking into consideration the elements of crossover point vectors, the parent chromosome sections between two consecutive crossover points are exchanged, thereby creating the offspring chromosomes. Figure 5 illustrates an example represented exclusively to describe the dynamics of the implemented code, but whose numbers do not match the specific problem dealt with in this study.

Journal Pre-proof

Figure 5: Representation of implemented binary x-point crossover A single generation of genetic algorithm involves the creation of a new population by executing selection and crossover several times. At each iteration, the selected parents and their offspring chromosomes are stored in a different array and finally sorted by applying the non-dominated and crowding distance principles. From the latter elitist ranking of chromosomes, the new population is created by selecting the first Pop_size solutions. Before inserting the chromosomes on the new population array, it is necessary to check if the offspring individuals meet the problem constraints. If the equations are not satisfied, the crossover is repeated until the offspring chromosomes are feasible. 5.4 Mutation To maintain the genetic diversity of chromosomes from one generation to the next one, a mutation operator is used. In the literature, several studies showed that mutation is an effective way for avoiding early convergence of an algorithm to a local optimal solution [94]. In this study, a counter is set up to define when the mutation has to be performed. In particular, whenever the counter is equal to the IndexMutation value, the mutation is performed and the counter is reset immediately. Therefore, chromosome mutation is performed at the steps (in terms of number of generations) of IndexMutation. When the mutation occurs, a comparison between a chromosome and all of the others is performed. If 40% of the genes are equal in the two chromosomes, the mutation is applied to the first chromosome. In particular, the genes involved in the mutation are found by creating the Mutation points row vector (1, q) consisting of integer random values in the range [1, genome length]. The array dimension q is equal to θ/2, i.e. the authors decide to mutate 50% of the individual genes. Likewise, the crossover operator, even for the mutation, is necessary for checking if the mutated chromosome is coherent with the problem constraints. If all of the equation constraints hold, the algorithm goes on comparing the next chromosome with the remaining ones; otherwise, the mutation

Journal Pre-proof is repeated on the same individual by mutating different genes (i.e. a new random mutation point vector is generated). Figure 6 shows a representative example of how the mutation is performed.

Figure 6: Representation of implemented binary mutation The method’s logical steps described above are shown in Figure 7 in which the flowchart of the implemented non-dominated sorting genetic algorithm is differentiated from the previous phase of the demand simulation. The flowchart includes also three check gateways, two of which refer to the feasibility check of the new solutions generated by the algorithm for the crossover and mutation operators, whereas the third ensures that the same parents are not selected.

Figure 7: Flowchart of simulation customer demand and genetic algorithm code

Journal Pre-proof

6. Solution selection with ELECTRE method ELECTRE is a family of MCDMs used for choosing, ranking, or sorting the best actions from a given set of alternatives according to several conflicting criteria [95; 96]. In this work, the ELECTRE III method is applied to find the ranking of solutions belonging to each Pareto front. MCDMs allow obtaining trade-off solutions with respect to several conflicting criteria. By means of the definition of proper thresholds, the ELECTRE III algorithm aims at determining first the outranking relations belonging to the [0,1] interval, according to which an alternative a outranks an alternative b if a is as good as b and good reasons to reject the last assertion do not exist. The computation of outranking relations is based on the conditions of concordance and discordance, which express the degree of preference and non-preference of a compared to b, respectively. The conversion of input data related to the decision maker evaluations into objective indicators requires the introduction of three thresholds for each criterion: indifference I, preference S, and veto V. The indifference threshold states the limit value of the difference between the scores of alternatives s(a) and s(b) on the same criterion below, in which the involved alternatives preserve the indifference. The preference threshold states the limit value of the difference between s(a) and s(b) above which a is preferred over b. The veto threshold states the limit value of the difference between s(a) and s(b) beyond which the gap of these scores is greatly excessive to be incompatible with the assertion ‘alternative a is strongly preferred to b’. For each criterion, the relation Ij ≤ Sj ≤ Vj holds. Let J = 1, 2,…, j…, G be the set of problem criteria, a and b the generic alternatives to be compared and ranked, and s j(a), s j(b) the score of alternatives a and b on criterion j, respectively. The necessary steps to implement the ELECTRE III algorithm are the following. The methodology starts with the computation of the marginal concordance indices cj (a, b) Ɐ j∈ J (Eq. (37)):

{

𝑐𝑗(𝑎,𝑏) = 1 𝑐𝑗(𝑎,𝑏) = 0 𝑐𝑗(𝑎,𝑏) =

𝑖𝑓 𝑠𝑗(𝑎) ≥ 𝑠𝑗(𝑏) ∧ 𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) ≤ 𝐼𝑗 𝑖𝑓 𝑠𝑗(𝑎) < 𝑠𝑗(𝑏) ∧ 𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) > 𝑆𝑗

𝑠𝑗(𝑎) ― 𝑠𝑗(𝑏) + 𝑆𝑗 𝑆𝑗 ― 𝐼𝑗

(37)

𝑖𝑓 𝑠𝑗(𝑎) < 𝑠𝑗(𝑏) ∧ 𝐼𝑗 < 𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) ≤ 𝑆𝑗

The marginal concordance matrices related to each criterion are then aggregated in a single comprehensive concordance index c*(a, b) (Eq. (38)), where 𝑤𝑗 is the weight-importance coefficient related to criterion j: 𝑐 ∗ (𝑎,𝑏) = ∑∀𝑗 ∈ 𝐽𝑤𝑗 ∙ 𝑐𝑗(𝑎,𝑏)

(38)

Journal Pre-proof The marginal discordance indices dj (a, b) Ɐ j∈ J are calculated by Eq. (39):

{

𝑑𝑗(𝑎,𝑏) = 0 𝑑𝑗(𝑎,𝑏) = 1 𝑑𝑗(𝑎,𝑏) =

𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) ― 𝑆𝑗 𝑉 𝑗 ― 𝑆𝑗

𝑖𝑓 𝑐𝑗(𝑎,𝑏) ≠ 0 𝑖𝑓 𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) ≥ 𝑉𝑗

(39)

𝑖𝑓 𝑠𝑗(𝑎) < 𝑠𝑗(𝑏) ∧ 𝑆𝑗 ≤ 𝑠𝑗(𝑏) ― 𝑠𝑗(𝑎) < 𝑉𝑗

Eq. (40) indicates the construction of the credibility index of the outranking relation σ(a, b) matrix:

{

𝜎(𝑎,𝑏) = 𝑐 ∗ (𝑎,𝑏) 𝑖𝑓 𝑑𝑗(𝑎,𝑏) = 0 ∀𝑗 ∈ 𝐽 ∗ 𝜎(𝑎,𝑏) = 𝑐 (𝑎,𝑏) 𝑖𝑓 ∃𝑗:𝑑𝑗(𝑎,𝑏) > 0 ∧ ∀𝑗 ∈ 𝐽 𝑑𝑗(𝑎,𝑏) < 𝑐 ∗ (𝑎,𝑏) 𝜎(𝑎,𝑏) = 𝑐 ∗ (𝑎,𝑏) ∙ ∏𝑗 ∈ 𝐽′

1 ― 𝑑𝑗(𝑎,𝑏) 1―𝑐

∗ (𝑎,𝑏)

𝑖𝑓 ∀𝑗 ∈ 𝐽′: 𝑑𝑗(𝑎,𝑏) ≥ 𝑐 ∗ (𝑎,𝑏)

(40)

where J’ denotes the set of criteria whose marginal discordance index dj is greater than or equal to the comprehensive concordance index c*(a, b). The computation of outranking relation indices σ(a, b) concludes the first phase of the methodology. The second one consists of implementing the ranking procedure by means of descending and ascending distillations. First, another threshold λ given by Eq. (41) must be set, which corresponds to the maximum value of σ(a, b) in the credibility matrix and from which a cut-off level of outranking λ’ is computed. The latter in Eq. (43) is equal to the difference between λ and a discrimination threshold that varies linearly with λ [97]. (41)

𝜆 = max 𝜎(𝑎,𝑏) ∀𝑎,𝑏

𝜑(𝜆) = 𝛼(𝜆) +𝛽 𝜆′ = 𝜆 ―

𝜑(𝜆)

(42) (43)

The cut-off level presented in Eq. (42) represents the minimum value at which the decision maker assigns a significant degree of credibility. To apply the distillation procedures, a Boolean matrix T(a, b) (Eq. (44)) is found: it contains a score of +1 for alternative a every time a outranks b; otherwise, a score equal to 0 is assigned.

{𝑇(𝑎,𝑏) = 1

𝑖𝑓 𝜎(𝑎,𝑏) ≥ 𝜆′ ∧ 𝜎(𝑎,𝑏) ― 𝜎(𝑏,𝑎) > 𝜑(𝜆) 𝑇(𝑎,𝑏) = 0 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

(44)

The final qualification score Q for an alternative is defined by the difference between the number of solutions it outranks and the number of solutions it is outranked by. The distillation is an iterative

Journal Pre-proof procedure that involves updating the Boolean matrix T and the subsequent qualification score computation and continues until all solutions are ranked. In particular, by considering the descending distillation, the alternative with the highest qualification score Q is assigned to the highest available rank and removed from the procedure by repeating this process for all remaining alternatives. On the other hand, the ascending distillation works exactly in the same manner; however, it starts from the alternative with the lowest qualification score Q by assigning it to the highest available rank. The final solution ranking is the result of the combination of the two distillation procedures. 7. Numerical experiment and computational results This section describes the numerical application of the 3E-FGVRP model discussed above and shows how the heuristic optimisation works in determining the Pareto fronts. The description of the entire model, implementation of NSGA-II, and relative MATLAB® analyses are presented. Table 4 provides all parameter values used in this study for executing the 3E-FGVRP optimisation. Table 4: Parameter values by 3E-FGVRP and NSGA-II Parameter

Value

A

42

C

9

D

4

E

56

Ee

0.03

Ef

0.042

F

7

FCk

5000

Fcn

6500

H

72

IndexMutation

20

L

72

MaxGeneration

5000

Mk

150

Pop_size

50

Qn

135

T

5

TCz

0.35

U

12

V

10

Journal Pre-proof α

0.8

β

0.6

γ

0.8

θ

13,110

Figure 8 shows the final Pareto fronts of the 3E-FGVRP of T periods. It can be noted that the solutions belonging to the represented Pareto fronts are characterised by very different values in terms of both costs and carbon emissions. These variances are due to the variety of customers’ demand input data in different periods, which resulted in different travelled arcs and routes in the 3E-FGVRP.

Figure 8: Pareto fronts of 3E-FGVRP in different time periods Finally, to choose a single solution of the supply chain for each simulated period, the application of the ELECTRE III multi-criteria method is proposed. The weights of the two objectives, i.e. total cost and carbon emission, are set as 0.4 and 0.6, respectively, because, as previously mentioned, the environmental impact of the supply chain is considered as a primary issue for its design. For both criteria, the thresholds required for implementing the ELECTRE III algorithm are defined as percentage values of differences between the scores of two alternatives (i.e. Pareto solutions). In particular, the indifference, preference and veto thresholds are set equal to 9%, 35% and 60% respectively for cost and carbon emission criteria. With relation to every considered scenario, the best solution is obtained by applying the two described distillation procedures. Obtained results are listed in Table 5. This table also indicates the number of

Journal Pre-proof vehicles and distribution centres involved in each supply chain. It is worth noting that for all scenarios, the number of vehicles is the same; moreover only in the scenario t5 all distribution centres are used. Table 5: Best solutions of Pareto fronts by means of ELECTRE III method Best solutions Carbon

Pop1

Total cost

t1

63,559

5985

7

3

t2

64,053

5904

7

3

t3

60,817

6860

7

3

t4

63,377

5058

7

3

t5

68,510

2777

7

4

emission

Vehicles Depots

To make the study more accurate and complete, a sensitivity analysis is finally conducted by repeating the heuristic optimisation procedure on three different random initial populations. For simplicity, only the Pareto fronts of period t2 are shown in Figure 9.

Figure 9: Pareto fronts of 3E-FGVRP in t2 for different initial populations The ELECTRE III algorithm, with the same thresholds and weights of the previous case, is also applied to the Pareto fronts shown in Figure 9 to determine the robustness degree of the NSGA-II solutions. Table 6 presents the objective function values of the best solutions related to the Pareto fronts of Figure 9. It is calculated that these solutions differ on average by 5% on the cost criteria; in contrast, as for the carbon emissions, they differ on average by 17%. Regarding the number of

Journal Pre-proof vehicles and utilised depots, the solutions lead to slightly different supply chains, even if two-thirds of the solutions converge to 7 vehicles and 4 depots. Therefore, it is possible to consider the latter as the definitive configuration for the period considered. The authors believe that these differences are relatively acceptable considering that they are found by means of a heuristic algorithm based on random mechanisms. Therefore, it can be concluded that the model is robust. Table 6: Best solutions of Pareto fronts obtained by three different populations Best solutions t2

Total cost

Carbon emission

Vehicles Depots

Pop1

64,053

5904

7

3

Pop2

60,740

6965

6

4

Pop3

66,118

5279

7

4

8. Conclusions In this study, the design of a three-echelon supply chain in a GVRP setting is performed. This involves the analysis of implementing a distribution chain related to ‘farm to table’ food products in a delimitated geographical area. This problem has been studied to date without taking into consideration the fuzzy customer demands. To make the sizing more accurate, the input data corresponding to the amount of food products requested by the customers are considered as triangular fuzzy numbers. Therefore, the resulting supply chain includes the uncertainty associated with the customers’ purchasing decision variability. The mathematical model of the 3E-FGVRP is developed by considering two objective functions, i.e. the total costs involved and the amount of CO2 emissions of the routes in each period. The fuzzy chance features are represented by means of problem constraints that are formulated and linearised by resorting to the credibility theory of fuzzy sets. Validating the model by the application of a case study and considering real values of the parameters used make this study an effective and truthful methodology for the design of a distribution chain. Moreover, this study proposes the integration of the green issue in the design phase, which makes the resulting chain consistent with current sustainable development and distribution themes. The optimisation is finally conducted by implementing NSGA-II in the MATLAB® environment, whose algorithm code is customised with respect to the conventional one. After having found the Pareto fronts, the MCDM ELECTRE III method is applied for determining the best solutions belonging to each of them. Finally, a sensitivity analysis is conducted by considering three different initial random populations for the heuristic optimisation procedure. It is found that the difference between the best solutions is acceptable; therefore, the methodology is considered robust.

Journal Pre-proof

Appendix Matrices of distances f1

f2

f3

f4

f5

f6

f7

d1

f1

d2

d3

d4

77

354 241 185

f2

62

68

310 184 128

f3

343 299

288

64

f4

279 222 165

350 144

f5

268 224

f6

287 174 104 113 210

233 115

70

f7

366 322

311

129 153

c1

c2

73 76

c3

236

42

168

217 169 277 187

156 161 127

c4

137 130

c5

c6

c7

c8

c1

c9

18

d1

d2

52

d3

d4

73 230 232 140

c2

224

296 86 193 217

c3

107 330

59 335 280 224

c4

211 96 316

280 17

c5

127 262 175 167

134 185 109 50

c6

140 182 246 88

c7

134 211 240 116 54

c8

261 162 366 67 216 137 164

331 85 115 172

c9

254 200 302 105 131 123 131 89

261 123 25 151

d1

d2

d3

86

99 122

210 106 101 42 36

204 132 107 10

d4

d1 d2

308

d3

238 117

d4

182 141 110

References [1]

Dantzing, G., Ramser, J. (1959). The truck dispatching problem. Management Science, 10(6), 80-91.

[2]

Orloff, C.S. (1974) A fundamental problem in vehicle routing. Networks, 4(1), 35-64.

[3]

Christofides, Nicos. (1976). Vehicle routing problem. Rev Fr Autom Inf Rech Oper, 10 (2), 5570.

[4]

Golden, B.L., Magnanti, T.L., Nguyen, H.Q. (1977). Implementing vehicle routing algorithms. Networks, 7(2), 113-148.

Journal Pre-proof [5]

Malandraki, Chryssi, Daskin, Mark S. (1992). Time dependent vehicle routing problems: Formulations, properties and heuristic algorithms. Transportation Science, 26(3), 185-200.

[6]

Solomon, Marius M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations Research, 35(2), 254-265.

[7]

Savelsbergh, M.W.P. (1985). Local search in routing problems with time windows. Annals of Operations Research, 4(1), 285-305.

[8]

Vidal, T., Crainic, T.G., Gendreau, M., Prins, C. (2013). A hybrid genetic algorithm with adaptive diversity management for a large class of vehicle routing problems with timewindows. Computers and Operations Research, 40(1), 475-489.

[9]

Teodorović, D., Pavković, G. (1996). The fuzzy set theory approach to the vehicle routing problem when demand at nodes is uncertain. Fuzzy Sets and Systems, 82(3), 307-317.

[10] Zheng, Y., Liu, B. (2006). Fuzzy vehicle routing model with credibility measure and its hybrid intelligent algorithm. Applied Mathematics and Computation, 176(2), 673-683. [11] Erbao, C., Mingyong, L. (2009). A hybrid differential evolution algorithm to vehicle routing problem with fuzzy demands. Journal of Computational and Applied Mathematics, 231(1), 302310. [12] Cao, E., Lai, M. (2010). The open vehicle routing problem with fuzzy demands. Expert Systems with Applications, 37(3), 2405-2411. [13] Brito, J., Martínez, F.J., Moreno, J.A., Verdegay, J.L. (2015). An ACO hybrid metaheuristic for close-open vehicle routing problems with time windows and fuzzy constraints. Applied Soft Computing Journal, 32, 154-163. [14] Shi, Y., Boudouh, T., Grunder, O. (2017) A hybrid genetic algorithm for a home health care routing problem with time window and fuzzy demand. Expert Systems with Applications, 72, 160-176. [15] Lyon, T.P., Maxwell, J.W. (2008). Corporate social responsibility and the environment: A theoretical perspective. Review of Environmental Economics and Policy, 2(2), 240-260. [16] Arminen, H., Puumalainen, K., Pätäri, S., Fellnhofer, K. (2018). Corporate social performance: Inter-industry and international differences. Journal of Cleaner Production, 177, 426-437. [17] De Lange, D.E. (2017). Start-up sustainability: An insurmountable cost or a life-giving investment? Journal of Cleaner Production, 156, 838-854. [18] Tsao, Y.C. (2015). Design of a carbon-efficient supply-chain network under trade credits. International Journal of Systems Science: Operations and Logistics. 2(3), 177-186.

Journal Pre-proof [19] Awasthi, A., Omrani, H. (2019). A goal-oriented approach based on fuzzy axiomatic design for sustainable mobility project selection. International Journal of Systems Science: Operations and Logistics, 6(1), 86-98. [20] Sayyadi, R., Awasthi, A. (2018). An integrated approach based on system dynamics and ANP for evaluating sustainable transportation policies. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2018.1554168. [21] Hao, Y., Helo, P., Shamsuzzoha, A. (2018). Virtual factory system design and implementation: integrated sustainable manufacturing. International Journal of Systems Science: Operations and Logistics, 5(2), 116-132. [22] Gharaei, A., Karimi, M., Hoseini Shekarabi, S.A. (2019). An integrated multi-product, multibuyer supply chain under penalty, green, and quality control polices and a vendor managed inventory with consignment stock agreement: The outer approximation with equality relaxation and augmented penalty algorithm. Applied Mathematical Modelling, 69, 223-254. [23] Sayyadi, R., Awasthi, A. (2018). A simulation-based optimisation approach for identifying key determinants for sustainable transportation planning. International Journal of Systems Science: Operations and Logistics, 5(2), 161-174. [24] Kazemi, N., Abdul-Rashid, S.H., Ghazilla, R.A.R., Shekarian, E., Zanoni, S. (2018). Economic order quantity models for items with imperfect quality and emission considerations. International Journal of Systems Science: Operations and Logistics, 5(2), 99-115. [25] Dubey, R., Gunasekaran, A., Sushil, Singh, T. (2015). Building theory of sustainable manufacturing using total interpretive structural modelling. International Journal of Systems Science: Operations and Logistics, 2(4), 231-247. [26] Lin, C., Choy, K.L., Ho, G.T.S., Chung, S.H., Lam, H.Y. (2014). Survey of Green Vehicle Routing Problem: Past and future trends. Expert Systems with Applications, 41(4), 1118-1138. [27] Jemai, J., Zekri, M., Mellouli, K. (2012). An NSGA-II algorithm for the green vehicle routing problem. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 7245 LNCS, 37-48. [28] Tiwari, A., Chang, P.-C. (2015). A block recombination approach to solve green vehicle routing problem. International Journal of Production Economics, 164, 379-387. [29] Koç, Ç., Karaoglan, I. (2016). The green vehicle routing problem: A heuristic based exact solution approach. Applied Soft Computing Journal, 39, 154-164. [30] Erdoĝan, S., Miller-Hooks, E. (2012). A Green Vehicle Routing Problem. Transportation Research Part E: Logistics and Transportation Review, 48(1), 100-114.

Journal Pre-proof [31] Bektaş, T., Laporte, G. (2011). The Pollution-Routing Problem. Transportation Research Part B: Methodological, 45(8), 1232-1250. [32] Demir, E., Bektaş, T., Laporte, G. (2012). An adaptive large neighborhood search heuristic for the Pollution-Routing Problem. European Journal of Operational Research, 223(2), 346-359. [33] Kramer, R.E., Subramanian, A., Vidal, T., Cabral, L.D.A.F. (2015). A metaheuristic approach for the Pollution-Routing Problem. European Journal of Operational Research, 243(2), 523539. [34] Hosseini-Nasab, H., Lotfalian, P. (2017). Green routing for trucking systems with classification of path types. Journal of Cleaner Production, 146, 228-233. [35] Micale, R., Marannano, G., Giallanza, A., Miglietta, P.P., Agnusei, G.P., La Scalia, G., (2019), Sustainable vehicle routing based on firefly algorithm and TOPSIS methodology, Sustainable Futures, 1 (2019), 100001. [36] Ting, C.-K., Liao, X.-L. (2013). The selective pickup and delivery problem: Formulation and a memetic algorithm. International Journal of Production Economics, 141(1), 199-211. [37] Gribkovskaia, I., Laporte, G., Shyshou, A. (2008). The single vehicle routing problem with deliveries and selective pickups. Computers and Operations Research, 35(9), 2908-2924. [38] Gutiérrez-Jarpa, G., Desaulniers, G., Laporte, G., Marianov, V. (2010). A branch-and-price algorithm for the Vehicle Routing Problem with Deliveries, Selective Pickups and Time Windows. European Journal of Operational Research, 2016(2), 341-349. [39] Aras, N., Aksen, D., Tuǧrul Tekin, M. (2011). Selective multi-depot vehicle routing problem with pricing. Transportation Research Part C: Emerging Technologies, 19(5), 866-884. [40] Soleimani, H., Chaharlang, Y., Ghaderi, H. (2018). Collection and distribution of returnedremanufactured products in a vehicle routing problem with pickup and delivery considering sustainable and green criteria. Journal of Cleaner Production, 172, 960-970. [41] Rabbani, M., Foroozesh, N., Mousavi, S.M., Farrokhi-Asl, H. (2019). Sustainable supplier selection by a new decision model based on interval-valued fuzzy sets and possibilistic statistical reference point systems under uncertainty. International Journal of Systems Science: Operations and Logistics, 6(2), 162-178. [42] Hoseini Shekarabi, S.A., Gharaei, A., Karimi, M. (2019). Modelling and optimal lot-sizing of integrated multi-level multi-wholesaler supply chains under the shortage and limited warehouse space: generalised outer approximation. International Journal of Systems Science: Operations and Logistics, 6(3), 237-257. [43] Rabbani, M., Hosseini-Mokhallesun, S.A.A., Ordibazar, A.H., Farrokhi-Asl, H. (2018). A hybrid robust possibilistic approach for a sustainable supply chain location-allocation network

Journal Pre-proof design. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2018.1506061. [44] Sarkar, S., Giri, B.C. (2018). Stochastic supply chain model with imperfect production and controllable defective rate. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2018.1536231. [45] Giri, B.C., Bardhan, S. (2014). Coordinating a supply chain with backup supplier through buyback contract under supply disruption and uncertain demand. International Journal of Systems Science: Operations and Logistics, 1(4), 193-204. [46] Hannan, M.A., Akhtar, M., Begum, R.A., Basri, H., Hussain, A., Scavino, E. (2018). Capacitated vehicle-routing problem model for scheduled solid waste collection and route optimization using PSO algorithm. Waste Management, 71, 31-41. [47] Gilardino, A., Rojas, J., Mattos, H., Larrea-Gallegos, G., Vázquez-Rowe, I. (2017). Combining operational research and Life Cycle Assessment to optimize municipal solid waste collection in a district in Lima (Peru). Journal of Cleaner Production, 156, 589-603. [48] Samanlioglu, F. (2013). A multi-objective mathematical model for the industrial hazardous waste location-routing problem. European Journal of Operational Research, 226(2), 332-340. [49] Faccio, M., Persona, A., Zanin, G. (2011). Waste collection multi objective model with real time traceability data. Waste Management, 31(12), 2391-2405. [50] Reed, M., Yiannakou, A., Evering, R. (2014). An ant colony algorithm for the multicompartment vehicle routing problem. Applied Soft Computing Journal, 15, 169-176. [51] Gharaei, A., Hoseini Shekarabi, S.A., Karimi, M. (2019). Modelling And optimal lot-sizing of the replenishments in constrained, multi-product and bi-objective EPQ models with defective products: Generalised Cross Decomposition. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2019.1574364. [52] Gharaei, A., Karimi, M., Hoseini Shekarabi, S.A. (2019). Joint Economic Lot-sizing in Multiproduct Multi-level Integrated Supply Chains: Generalized Benders Decomposition. International

Journal

of

Systems

Science:

Operations

and

Logistics,

DOI:

10.1080/23302674.2019.1585595. [53] Gharaei, A., Hoseini Shekarabi, S.A., Karimi, M., Pourjavad, E., Amjadian, A. (2019). An integrated stochastic EPQ model under quality and green policies: generalised cross decomposition under the separability approach. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2019.1656296. [54] Sasikumar, P., Kannan, G. (2008). Issues in reverse supply chains, part II: Reverse distribution issues - an overview. International Journal of Sustainable Engineering, 1(4), 234-249.

Journal Pre-proof [55] Lin, C., Choy, K.L., Ho, G.T.S., Ng, T.W. (2014). A Genetic Algorithm-based optimization model for supporting green transportation operations. Expert Systems with Applications, 41(7), 3284-3296. [56] Kim, H., Yang, J., Lee, K.-D. (2011). Reverse logistics using a multi-depot vrp approach for recycling end-of-life consumer electronic products in South Korea. International Journal of Sustainable Transportation, 5(5), 289-318. [57] Habibi, M.K.K., Battaïa, O., Cung, V.-D., Dolgui, A. (2017). Collection-disassembly problem in reverse supply chain. International Journal of Production Economics, 183, 334-344. [58] Habibi, M.K.K., Battaïa, O., Cung, V.-D., Dolgui, A. (2017). An efficient two-phase iterative heuristic for Collection-Disassembly problem. Computers and Industrial Engineering, 110, 505-514. [59] Lin, C., Choy, K.L., Ho, G.T.S., Ng, T.W. (2014). A Genetic Algorithm-based optimization model for supporting green transportation operations. Expert Systems with Applications, 41(7), 3284-3296. [60] Dell’Amico, M., Righini, G., Salani, M. (2006). A branch-and-price approach to the vehicle routing problem with simultaneous distribution and collection. Transportation Science, 40(2), 235-247. [61] Tasan, A.S., Gen, M.(2012). Genetic algorithm based approach to vehicle routing problem with simultaneous pick-up and deliveries. Computers and Industrial Engineering, 62(3), 755-761. [62] Xiao, Y., Konak, A. (2018). A genetic algorithm with exact dynamic programming for the green vehicle routing & scheduling problem. Journal of Cleaner Production, 167, 1450-1463. [63] Qiu, Y., Wang, L., Fang, X., Pardalos, P.M., Goldengorin, B. (2018). Formulations and branchand-cut algorithms for production routing problems with time windows. Transportmetrica A: Transport Science, 1-22. [64] Feillet, D., Dejax, P., Gendreau, M., Gueguen, C. (2004). An exact algorithm for the elementary shortest path problem with resource constraints: Application to some vehicle routing problems. Networks, 44(3), 216-229. [65] Laporte, G., Louveaux, F.V., Van Hamme, L. (2002). An integer L-shaped algorithm for the capacitated vehicle routing problem with stochastic demands. Operations Research, 50(3), 415423. [66] Baldacci, R., Hadjiconstantinou, E., Mingozzi, A. (2004). An exact algorithm for the capacitated vehicle routing problem based on a two-commodity network flow formulation. Operations Research, 52(5), 723-738.

Journal Pre-proof [67] Valério De Carvalho, J.M. (2002). LP models for bin packing and cutting stock problems. European Journal of Operational Research, 141(2), 253-273. [68] Tarantilis, C.D., Kiranoudis, C.T. (2001). A meta-heuristic algorithm for the efficient distribution of perishable foods. Journal of Food Engineering, 50(1), 1-9. [69] Baker, B.M. Ayechew, M.A. (2003). A genetic algorithm for the vehicle routing problem. Computers and Operations Research, 30(5), 787-800. [70] Hiassat, A., Diabat, A., Rahwan, I. (2017). A genetic algorithm approach for locationinventory-routing problem with perishable products. Journal of Manufacturing Systems, 42, 93103. [71] Koç, Ç., Bektaş, T., Jabali, O., Laporte, G. (2016). The fleet size and mix location-routing problem with time windows: Formulations and a heuristic algorithm. European Journal of Operational Research, 248(1), 33-51. [72] Karakatič, S., Podgorelec, V. (2015). A survey of genetic algorithms for solving multi depot vehicle routing problem. Applied Soft Computing Journal, 27, 519-532. [73] Liu, R., Xie, X., Augusto, V., Rodriguez, C. (2013). Heuristic algorithms for a vehicle routing problem with simultaneous delivery and pickup and time windows in home health care. European Journal of Operational Research, 230(3), 475-486. [74] Vidal, T., Crainic, T.G., Gendreau, M., Lahrichi, N., Rei, W. (2012). A hybrid genetic algorithm for multidepot and periodic vehicle routing problems. Operations Research, 60(3), 611-624. [75] Ghoseiri, K., Ghannadpour, S.F. Multi-objective vehicle routing problem with time windows using goal programming and genetic algorithm. Applied Soft Computing Journal, 10(4), 10961107. [76] Shah, N.H., Chaudhari, U., Cárdenas-Barrón, L.E. (2018). Integrating credit and replenishment policies for deteriorating items under quadratic demand in a three echelon supply chain. International

Journal

of

Systems

Science:

Operations

and

Logistics,

DOI:

10.1080/23302674.2018.1487606. [77] Giri, B.C., Masanta, M. (2018). Developing a closed-loop supply chain model with price and quality dependent demand and learning in production in a stochastic environment. International Journal of Systems Science: Operations and Logistics, DOI: 10.1080/23302674.2018.1542042. [78] Yin, S., Nishi, T., Zhang, G. (2016). A game theoretic model for coordination of single manufacturer and multiple suppliers with quality variations under uncertain demands. International Journal of Systems Science: Operations and Logistics, 3(2), 79-91.

Journal Pre-proof [79] Duan, C., Deng, C., Gharaei, A., Wu, J., Wang, B. (2018). Selective maintenance scheduling under stochastic maintenance quality with multiple maintenance actions. International Journal of Production Research, 56(23), 7160-7178. [80] Zadeh, L.A. (1965). Fuzzy sets. Information and Control, 8(3), 338-353. [81] Zadeh LA. (1978). Fuzzy sets as a basis for a theory of possibility. Fuzzy Sets Syst, 1(1), 3-28. [82] Liu, Baoding (2004). Uncertainty Theory: An Introduction to its Axiomatic Foundations. [83] Liu, B., Liu, Y.-K. (2002). Expected value of fuzzy variable and fuzzy expected value models. IEEE Transactions on Fuzzy Systems, 10(4), 445-450. [84] Dubois, D., Prade, H. (1978). Operations on fuzzy numbers. International Journal of Systems Science, 9(6), 613-626. [85] www.aci.it [86] Adiba, E.B.E.I., Elhassania, M., Ahemd, E.H.A. (2013). A Hybrid Ant Colony System for Green capacitated vehicle routing problem in sustainable transport. Journal of Theoretical and Applied Information Technology, 54(2), 198-208. [87] Pan, S., Ballot, E., Fontane, F. (2013). The reduction of greenhouse gas emissions from freight transport by pooling supply chains. International Journal of Production Economics, 143(1), 8694. [88] Palmer A. (2007). The development of an integrated routing and carbon dioxide emissions model for goods vehicles. School of Management Cranfield University. [89] Deb K, Agrawal S, Pratap A, Meyarivan T. A fast elitist non-dominated sorting genetic algorithm for multi-objective optimization: NSGA-II. Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics) 2000; 1917: 849-858. [90] Zitzler, E., Deb, K., Thiele, L. (2000). Comparison of multiobjective evolutionary algorithms: empirical results. Evolutionary computation, 8(2), 173-195. [91] Deb, K., Pratap, A., Agarwal, S., Meyarivan, T. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation, 6(2), 182-197. [92] Seshadri, Aravind. A FAST ELITIST MULTIOBJECTIVE GENETIC ALGORITHM: NSGAII. [93] Whitley, D. (1994). A genetic algorithm tutorial. Statistics and Computing, 4(2), 65-85. [94] Salhi, S., Gamal, M.D.H. (2003). A Genetic Algorithm Based Approach for the Uncapacitated Continuous Location-Allocation Problem. Annals of Operations Research, 123(1-4), 203-222. [95] Marzouk, M.M. (2011). ELECTRE III model for value engineering applications. Automation in Construction, 20(5), 596-600.

Journal Pre-proof [96] Miglietta, P.P., Micale, R., Sciortino, R., Caruso, T., Giallanza, A., La Scalia, G., (2019), “The sustainability of olive orchard planting management for different harvesting techniques: An integrated methodology”, Journal of Cleaner Production, 238, 117989. [97] Figueira, J., Mousseau, V., Roy, B. (2005). ELECTRE methods. International Series in Operations Research and Management Science, 78, 133-162.

Journal Pre-proof CREDIT AUTHOR STATEMENT Antonio Giallanza: Conceptualization, Methodology, Data curation, Writing-Original draft preparation, Visualization, Investigation, Supervision, Validation, Reviewing and Editing. Gabriella Li Puma: Conceptualization, Methodology, Data curation, Writing-Original draft preparation, Visualization, Investigation, Software.

Journal Pre-proof

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: