Alessandro Di Pretoro

Stéphane Negny

Ludovic Montastruc
email: ludovic.montastruc@ensiacet.fr

The Traveling Deliveryman Problem under Uncertainty: Fundamentals for Flexible Supply Chains

Keywords: uncertainty modelling, supply chain, exibility, traveling deliveryman problem, operational research

I

n recent years, supply chain optimization has become a major topic in operational research. Huge eorts were made to outline eective methodologies both from computational and nal solution quality perspective. However, besides being optimal, one major need of supply chains is exibility and adaptability to disturbances. This research work aims at the denition of a thorough procedure for exible supply chain design able to exploit deterministic and stochastic exibility indices. The Traveling Deliveryman Problem was then selected as case study since it represents the most basic operational research problem. The adaptation of exibility assessment methodologies, primarily conceived for unit operations, provides good results and allows to properly identify criticalities. In addition, it proves that the optimal solutions under uncertainty are conned to a limited subset among million possible alternatives according to the expected deviation nature. This procedure is then worth to be extended to more complex systems in future studies.

Highlights:

1. This study outlines a general procedure for supply chain exibility assessment that can be scaled up to more complex systems 2. The indicators conventionally used for unit operations have been adapted to the operational research providing good results 3. Detecting the most critical parameters allows to have more reliable cost expectations and to take dedicated measures

Introduction

During the last decades, exibility has seen a renewed attention as a fundamental requirement for the majority of production systems given the fast changing reality surrounding us. An important boost in this direction was given by the sustainability goals [1,[START_REF] Montastruc | Renewables 2020 global status report[END_REF] that forced the process industry to deal with green raw materials and renewable energy resources that show an intrinsically unstable nature.

In the Process Systems Engineering domain it was properly dened along with its related indicators at the beginning of 80's by Grossmann and Morari [1]. Since then, exibility analysis methodologies were mainly addressed to the industrial domain in general both to assess the exibility of an existing system (exibility index problem) or to design it for a given exibility requirement (exible design problem). A substantial number of studies about exible processes can be found in literature ranging from heat exchanger networks [START_REF] Swaney | An index for operational exibility in chemical process design. Part I: Formulation and theory[END_REF][START_REF] Saboo | Design of Resilient Processing Plants .8. a Resilience Index for Heat-Exchanger Networks[END_REF][START_REF] Pistikopoulos | A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters[END_REF][START_REF] Floudas | Global optimization in de-sign under uncertainty: feasibility test and exibility index problems[END_REF] to unit operations [START_REF] Hoch | Evaluation of design exibility in distillation-columns using rigorous models[END_REF][START_REF] Di Pretoro | Flexibility assessment of a biorenery distillation train: Optimal design under uncertain conditions[END_REF] or chemical reactors [START_REF] Huang | Operation feasibility analysis of chemical batch reaction systems with production quality consideration under uncertainty[END_REF]. Moreover, even if in a lower amount, the same indicators have been applied also to case studies concerning aspects other than equipment design such as maintenance scheduling [START_REF] Di Pretoro | Optimal cleaning cycle scheduling under uncertain conditions: A exibility analysis on heat exchanger fouling[END_REF],

process safety [START_REF] Castillo-Borja | A resilience index for process safety analysis[END_REF] or again economics [START_REF] Tan | A shock absorption index for inoperability input-output models[END_REF].

Although design under uncertainty and exibility were originally conceived for process related purposes, in the last decade, they are getting more and more attention from other domains as well. The operational research eld, logistics in particular, is indeed very sensitive to the impact of uncertainty on the quality of the related optimal results. That's why in recent years the uncertainty implications [START_REF] Sreedevi | Uncertainty and supply chain risk: The moderating role of supply chain exibility in risk mitigation[END_REF] and exibility indicators [START_REF] Behzadi | On metrics for supply chain resilience[END_REF] have become a major topic of research and innovation in this domain as well.

In fact, from a supply chain perspective, a more reliable assessment of the actual costs in case perturbations should occur as well as the identication of the proper actions to take in order to improve its exibility could be an important added value for the industrial domain from a logistics perspective. A supply chain is dened as a network within an organization or between multiple organizations that involves the procurement of raw materials, conversion from raw materials to nal products, and distribution of nal products to markets [START_REF] Simchi-Levi | Designing and managing the supply chain: concepts, strategies, and case studies[END_REF]. Thus, besides resources and products ows, it also involves economic and information aspects. As a consequence of this denition, supply chains are thus aected by every disturbance related to each of them.

In fact, in order to deal with a more and more volatile market demand, with variable resources to be exploited and with production changeover [START_REF] Simangunsong | Supply-chain uncertainty: a review and theoretical foundation for future research[END_REF][START_REF] Knight | Supply market uncertainty: Exploring consequences and responses within sustainability transitions[END_REF], accounting for uncertainty in supply chain management could be of particular interest in a short term perspective yet. To be ready soon to take up these challenges that are going to become even more complex during the next years, prompt and eective research works in this direction are required. In fact, this eld still has considerable room for improvement given the small amount of research works currently available in literature that couples exibility indices and supply chain management [START_REF] Govindan | Supply chain network design under uncertainty: A comprehensive review and future research directions[END_REF][START_REF] Cole | Supply Chain Optimization under Uncertainty. Supply Chain Design for Optimum Performance[END_REF][START_REF] Sahay | Flexibility assessment and risk management in supply chains[END_REF][START_REF] Baniasadi | A transformation technique for the clustered generalized traveling salesman problem with applications to logistics[END_REF]. However, an interesting application of exibility indices for supply chain management by Wang et al. [START_REF] Mastragostino | Flexibility analysis of process supply chain networks[END_REF] was proposed in 2016. In their work the authors carry out the exibility assessment of a chemicals supply network towards manufacturing sites including the distribution to customers accounting for demand, supply, yield and market exibility.

In this research we want to take a step back with respect to that article in order In particular, based on the analysis of Di Pretoro et al. (2021) [START_REF] Di Pretoro | Flexibility Analysis in Supply Chain Management: Application to the Traveling Salesman Problem[END_REF], a simple but appropriate case study will replace the complex chemicals supply chain one in order to ease the calculation aspects and emphasize the understanding of the uncertainty impact on the nal result. However, despite its lower computational eort, the operational research problem we are proposing carries with it all the characteristic features of more complex optimization examples and its choice does not restrain the validity of the obtained results.

In the following section the Traveling Deliveryman Problem case study will be presented in detail. In particular, the reasons behind this choice as well as the rigorous mathematical formulation of the optimization problem are thoroughly discussed.

2

Case study

The selected case study for this research work is a particular formulation of the Traveling Deliveryman Problem (hereafter TDP) [START_REF] Fischetti | The Delivery Man Problem and Cumulative Matroids[END_REF][START_REF] Méndez-Díaz | A new formulation for the Traveling Deliveryman Problem[END_REF] sometimes also mentioned as Traveling Repairman Problem or, with minor modications, as Vehicle Routing Problem (hereafter VRP).

This operational research problem belongs to the family of Minimum Cost

Hamiltonian Path Problem whose simplest formulation is given by the Traveling Salesman Problem (hereafter TSP) which aims at the search of the minimum total traveling distance between n cities as shown in Figure 1 for a better visualization. Besides the TDP, the TSP itself accounts for several similar formulations that mainly dier from each other because of the included constraints or some modications in the associated cost function without altering the nature of the problem itself. By adapting these modications to the TDP, the analogous versions of the TDP can be formulated as well, such as the TDP with time windows [START_REF] Tsitsiklis | Special cases of traveling salesman and repairman problems with time windows[END_REF][START_REF] Heilporn | The Delivery Man Problem with time windows[END_REF], multiple pickup and delivery [START_REF] Cheang | Multiple pickup and delivery traveling salesman problem with last-in-rst-out loading and distance constraints[END_REF], multiple objectives [START_REF] Wang | Uncertain multiobjective traveling salesman problem[END_REF][START_REF] Bock | Combining Traveling Salesman and Traveling Repairman Problems: A multi-objective approach based on multiple scenarios[END_REF] or stochastic travel times [START_REF] Laporte | The Vehicle Routing Problem with Stochastic Travel Times[END_REF][START_REF] Errico | A priori optimization with recourse for the vehicle routing problem with hard time windows and stochastic service times[END_REF].

The TSP represents the most known NP-hard problem in combinatorial optimization and, during the last century, its most eective solution algorithm has caught the attention of scientists from several research elds. Indeed, dierent procedures have been outlined for its exact solution and even more for its approximation or heuristic algorithms aimed at a substantial reduction of the involved computational eort when the size of the system, i.e. n, starts growing.

A detailed discussion about TSP solution algorithms can be found in Applegate et al. ( 2006) [START_REF] Applegate | The Traveling Salesman Problem: A Computational Study[END_REF] and in other literature works [START_REF] Roy | A Noble Genetic Algorithm to Solve a Solid Green Traveling Purchaser Problem with Uncertain Cost Parameters[END_REF][START_REF] Nechita | On the Resilience of Ant Algorithms. Experiment with Adapted MMAS on TSP[END_REF].

However, the computational eort is not the aspect this study focuses on.

The reason behind the TDP choice as a case study is that, although it is one of the simplest formulation in operational research problems, it shows all their features and critical aspects and, thus, it is a suitable example to set the basis for a exibility analysis approach of general validity in this research eld. In fact, although its low complexity, this problem allows to consider both locations and distances of dierent customers and to account for a cost function related to the load to be delivered along with their related uncertainties. Moreover, by adjusting these parameters more complex operational research problems can be derived without compromising the validity of the proposed procedure.

According to the purpose of this work then, the TDP formulation that will be addressed states as follows:

Given a set of n cities, a xed starting point representing the depot vertex and n loads to be respectively delivered to each of them whose weight represents a proportional contribution to the traveling time, the path minimizing the total traveling time should be found.

In order to avoid misunderstandings, a proper and cautious notation is required for the problem discussion that will follow in the next sections. As regards the TDP presented in this study, in order to make its mathematical formulation as intuitive as possible, cities are labeled both with an identifying letter, that is related to the aleatory point generation only, and with subscript numbers according to the traveling sequence (an example is given in Figure 1). For the unloaded traveling time between two cities i and j the expression t ij will be used while for the i-th city delivery l i represents the corresponding load. Every time a city is reached, the corresponding delivery is accomplished and, thus, the related load is discharged. Therefore, the load to be delivered to the j-th city after j iterations is given by:

l j = l tot - j-1 i=0 l i (2.1)
In order to have a comparable impact of loads and traveling distances on the optimal solution, the l vector is normalized so that l tot = n. Nonetheless, the implemented code presented in this work can be easily modied in order to change the l tot relative weight assumption.

The cost function to minimize resulting from this notation can be then writ-

C tot = n i=0,j̸ =i t ij • (1 + l j ) (2.2)
To carry out this study, one depot vertex and 10 cities have been used in order to ensure a complexity suciently high not to fall in particular cases and to keep a reasonable computational time. The depot vertex was xed at the origin of the axes and the ten cities were randomly located with coordinates in the 2D space within the range [-5, 5]. The load of the deliveries corresponding to each city have been randomly generated as well and normalized afterwards so that the initial total load l tot weigh as much as the unloaded traveling distance if the cities were uniformly aligned over the interval. To better understand on this case study hypothesis, from a mathematical point of view it can be formulated as:

n i=1 l i = n • ∆x (2.3)
where ∆x is the discretization step of the 2D space.

The details concerning the particular case study addressed in this research work are provided in Table 1. However, it is worth remarking once more that the validity of the procedure is not aected by the specic selection of cities location and delivery loads and the methodology can be applied to whatever number of cities and process parameters values.

Another important aspect that deserves particular attention for the purpose of this research work is the uncertainty characterization. This problem can be addressed by answering to two main questions, namely:

• what are the parameters of the system to be considered as uncertain;

• in case of stochastic exibility analysis, what is the probability distribution describing the likelihood that these parameters take a given value.

For the TDP case study, the main process parameters are the traveling distance between each pair of cities t ij = t ji and the load l i to be delivered to each of them. In this work both of them will be considered as uncertain parameters since they both play a key role in the operational research problem. However, for n cities, if all the distances and loads were assumed as uncertain, the uncertain domain would have a dimension equal to all the possible pairs plus the number of loads, i.e. n! 2!(n-2)! + n, that is a huge space with respect to the available computational power even for a limited value of n and, more importantly, that does not serve the purpose of the study. Therefore, in order to reduce the computational eort and, in particular, to emphasize the implication of each parameter on the nal optimum, the most critical of them from a exibility perspective have been considered as uncertain.

As better explained in the result section, the distance between the two farthest and the two closest cities along their connection direction, i.e. the longest and shortest travel paths, and the biggest and smallest loads are the system variables that show a more considerable impact on the best traveling path. In this research work, the two type of uncertainties are addressed separately in order to be able to show the obtained results on a visible geometrical space (3D or lower). However, the methodology keeps being valid even in case of simultaneous distance and load uncertainty and the exibility related charts can be nevertheless outlined according to the same procedure.

As regards the probability distribution associated to each of them for the stochastic exibility assessment further details are provided in the dedicated section that follows.

Methodology

This section is dedicated to the mathematical methodologies involved in this research work to solve the coupled Traveling Deliveryman Problem and its exibility problem.

The rst part refers to the fundamental properties of the exibility analysis and to the detailed description of the indicators that are employed in this paper to carry out the TDP results assessment under uncertain operating conditions.

In this section two main categories of exibility indices are presented along with the related challenges from a computational point of view.

On the other hand, the second section focuses on the algorithm used in this case study to solve the TDP problem as well as on the implications related to the selected uncertain variables.

Flexibility indices

During the last decades, exibility has become one of the topics of major concern in the Industrial Engineering domain and exible design is the main goal of process industry in order to cope with variable operating conditions related to the use of sustainable raw materials and energy sources or to the volatile market demand and prices. This property is dened as the capability of a system to accommodate a set of uncertain parameters without becoming infeasible.

In order to quantify the exibility of a given system with respect to uncertain parameter perturbations a suitable indicator is then required. For this purpose, several exibility indices have been proposed in literature since 1983 [START_REF] Swaney | An index for operational exibility in chemical process design. Part I: Formulation and theory[END_REF][START_REF] Saboo | Design of Resilient Processing Plants .8. a Resilience Index for Heat-Exchanger Networks[END_REF][START_REF] Pistikopoulos | A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters[END_REF][START_REF] Lai | Measurement of Plant Flexibility[END_REF].

As discussed by Di Pretoro et al. (2019) [START_REF] Di Pretoro | Flexibility analysis of a distillation column: Indexes comparison and economic assessment[END_REF], despite the amount of existing indicators, they can be classied into two main categories according to the way uncertainty is characterized, namely deterministic and stochastic indices.

On the one hand, the rst ones assess the exibility of a system as a measure of the maximum perturbation magnitude that can be withstood by the system without becoming infeasible. On the other hand, the stochastic indices are based on a probability function that associates each perturbed operating condition to its occurrence likelihood. Therefore, they evaluate the probability that the system is able to accommodate the eventual disturbances or, equivalently, the fraction of all the possible disturbances accommodated by the system.

In recent literature, a wide range of applications involving the various available indices can be found with a particular focus on industrial systems, while they have rarely been exploited for the description of other kind of systems. In this research work the most established index for each of the two categories will be used in order to show analogies and dierences related to the impact of the disturbance likelihood on the TDP nal results. The Swaney and Grossmann [START_REF] Swaney | An index for operational exibility in chemical process design. Part I: Formulation and theory[END_REF] deterministic exibility index F SG and the Stochastic Flexibility index SF

proposed by Pistikopoulos and Mazzuchi [START_REF] Pistikopoulos | A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters[END_REF] are then presented and discussed in deepth here below in order to provide a detailed understanding of their mathematical formulation and a clear view of the exibility analysis expected outcome.

Deterministic exibility

The rst exibility index and the related problem was introduced in 1983 by Swaney and Grossman [START_REF] Swaney | An index for operational exibility in chemical process design. Part I: Formulation and theory[END_REF] and has been widely used over the last decades to assess the exibility of various systems related to the process industry ranging from heat exchanger networks to unit operations [START_REF] Swaney | An index for operational exibility in chemical process design. Part I: Formulation and theory[END_REF][START_REF] Saboo | Design of Resilient Processing Plants .8. a Resilience Index for Heat-Exchanger Networks[END_REF][START_REF] Pistikopoulos | A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters[END_REF][START_REF] Floudas | Global optimization in de-sign under uncertainty: feasibility test and exibility index problems[END_REF][START_REF] Hoch | Evaluation of design exibility in distillation-columns using rigorous models[END_REF][START_REF] Di Pretoro | Flexibility assessment of a biorenery distillation train: Optimal design under uncertain conditions[END_REF][START_REF] Huang | Operation feasibility analysis of chemical batch reaction systems with production quality consideration under uncertainty[END_REF][START_REF] Di Pretoro | Optimal cleaning cycle scheduling under uncertain conditions: A exibility analysis on heat exchanger fouling[END_REF][START_REF] Castillo-Borja | A resilience index for process safety analysis[END_REF][START_REF] Lai | Measurement of Plant Flexibility[END_REF].

From an operational point of view, it is dened as the maximum fraction of the expected deviation of all the uncertain parameters at once that can be accommodated by the system.

The mathematical formulation of the Swaney and Grossman exibility index problem states as:

F SG = max δ (3.1) s.t. max θϵT (δ) min z max jϵJ f j (d, z, θ) ≤ 0 (3.2)
where:

• T (δ) is the hyper-rectangle described by the equation:

T (δ) = {θ : θ N -δ • ∆θ -< θ < θ N + δ • ∆θ + } (3.3)
• θ N represents the nominal vector of the uncertain parameters;

• ∆θ + , ∆θ -are the expected deviations in the positive and negative direction respectively;

• f j are the constraint functions that outline the feasible region;

• d and z are the design and the control (i.e. manipulated) variables respectively.

From a geometrical point of view, it corresponds to nding the scale factor for the largest hyperrectangle, whose sides are proportional to the i-th uncertain parameter expected deviation, that can be inscribed inside the feasible domain as shown by the blue rectangle in Figure 2. Furthermore, the corresponding deviation magnitude is represented in red.

According to the specic nature of the constraints outlining the feasible space in the uncertain domain, some strategies to reduce the complexity of the optimization problem solution can be implemented. In particular, for constraints

θ k = θ N + δ • ∆θ k (3.7)
Moreover, whether these hypotheses are satised, the solution lies at a vertex of the hyperrectangle. This means that the optimization problem solution could be substantially simplied as the evaluation of the system feasibility in correspondence of each hyperrectangle vertices for the corresponding δ, i.e. 2 n calculations. This procedure allows to considerably circumvent the computational eort required for the explicit solution of the min-max problem.

However, if this property of the feasible region cannot always be ensured, certain types of non-convex domains might lead to nonvertex solutions [START_REF] Grossmann | Active constraint strategy for exibility analysis in chemical processes[END_REF].

In order to have a more intuitive understanding of this behaviour let us refer once again to Figure 2. The result of the vertex analysis for the non-convex domain is shown by the dashed blue rectangle along with the corresponding deviation magnitude in red. As it can be immediately pointed out, in this case, the vertex analysis considerably overestimates the system exibility with respect to its actual value. Therefore, for a general set of constraints concerning a exibility index problem with n uncertain parameters, the feasibility over the entire external hypersurface of the n-dimensional rectangle at least should be assessed for each value of δ. Unfortunately, for the TDP it can be easily shown that, dierently from case studies involving process systems and equipment, the feasible space is more likely to be non-convex rather than convex. The reasons behind this behaviour will be explained in detail in the following results section. The main consequence of this phenomenon is then that, in order to carry out the exibility assessment related to this case study, the advantages of the vertex analysis can seldom be exploited.

Stochastic exibility

After a few years, Pistikopoulos and Mazzuchi (1990) [START_REF] Pistikopoulos | A Novel Flexibility Analysis Approach for Processes with Stochastic Parameters[END_REF] suggested to integrate the exibility analysis with the idea of perturbation probability in order to emphasize those disturbances that are more likely to occur and mitigate the impact of those who rarely take place on the calculated index. Their work led then to the conception of a new exibility indices category named as stochastic whose mathematical formulation states as follows.

Given the uncertain parameters Probability Distribution Function (hereafter PDF) P (θ) and the feasible region Ψ(d, z, θ) ≤ 0, the Stochastic Flexibility index SF can be dened as:

SF = ¢ Ψ P (θ) • dθ (3.8)
To better visualize its meaning, in Figure 3 Unlike the deterministic case, the calculation of the stochastic exibility indices does not substantially vary from the convex to the non-convex domain shape since, even in the rst case, data from the entire feasible space need to be used. This higher results accuracy is nevertheless paid not only by the higher computational eort but also by the higher amount of information needed, namely the probability characterization of all the points of the uncertain domain.

In this work, since no information is available about the perturbation occurrence, the stochastic exibility assessment was carried out with two well-known probability distribution functions. The rst one is the most general Gaussian (or Normal) multivariate PDF while the second one is a skewed distribution, namely the Beta PDF, shown respectively in Figures 4a and4b for the bidimensional domain. These two distributions have been selected to represent two common system behaviours: the rst one corresponds to the case where the disturbance likelihood is uniformly distributed with respect to the nominal operating conditions (symmetric PDF) while the second one to those cases when positive (or negative) deviations are more likely to occur. This impact of the PDF skewness on the resulting stochastic exibility index has been also pointed out in previous research works [START_REF] Di Pretoro | Optimal cleaning cycle scheduling under uncertain conditions: A exibility analysis on heat exchanger fouling[END_REF] and already proved to be something worth to account for.

Both the Gaussian and the Beta PDFs belong to the two parameter distribution family, that means they require two conditions to be uniquely dened [START_REF] Severini | Elements of Distribution Theory[END_REF]. The conditions used in this research work for the denition of the Gaussian PDF are:

• the mean value µ, that is also the mode value for the selected PDF, equal to the nominal operating conditions;

• the variance value σ selected in a way that the 99 % of the cumulative PDF lies within the expected deviation range.

For the Normal PDF the mean value and the variance are explicit in its analytical expression while for the Beta one they should be derived as a function of the so called shape (α) and scale (β) parameters. While a criterion analogous to the variance one was used for the scale parameter, the shape parameter was selected in order to shift the PDF mean value so that the 50 % of its CDF lies in one quarter of the uncertain domain. Although this can be seen as an arbitrary choice, it serves the purpose of the study to emphasize the dierence between skewed and symmetric uncertainty characterization. This means that, for PDF with higher skewness, the dierences highlighted in the results section between this Gaussian and Beta distributions, identied by means of the selected (µ, σ) and (α, β) pairs respectively, will be more evident.

Finally, it worth remarking that, even if the stochastic exibility assessment results will necessarily reect both qualitatively and quantitatively the properties of the selected PDF, the proposed procedure keeps being of general validity.

The TDP solution algorithm

Before presenting the detailed discussion of the TDP solution under uncertain operating conditions, the solution procedure is nevertheless worth a few remarks. As already explained in the case study section, the TDP is an NP-hard problem and exact algorithms, i.e. algorithms that always solve the optimization problem to optimality, can solve it in exponential time. The most trivial exact solution of this problem is obtained by enumeration, i.e. the minimum search among all possible permutations. This methodology nevertheless experiences a factorial growth of the computational time and it results impractical for a number of cities higher than 20, that is not the case of the presented case study.

Finally, a wide range of heuristic and approximation algorithms can be used [START_REF] Applegate | The Traveling Salesman Problem: A Computational Study[END_REF] to obtain very good solutions, 2-3 % error, in a quicker computational time and this is still an active research eld. Anyways, due to the limited number of cities in the presented case study, the best solution strategy to reduce the computational eort for the TDP solution under uncertain conditions will not be investigated since it does not serve the purpose of this research work, that is the detailed analysis of the exibility assessment results for dierent sets of uncertain parameters. However, it could be possibly addressed in later studies involving larger sets of points or more complex supply chain management case studies.

For the specic example used for this study, the solution of a single TDP under nominal operating conditions was performed by means of the enumeration algorithm implemented in a Matlab code and run on an Intel ® Core i7-7500 CPU with 2.7 GHz. The resulting computational time was in the order of 3 seconds.

Results

In this section the results of the exibility assessment carried out on the TDP case study are presented in the corresponding subsections classied according to the parameters that have been considered as uncertain. Due to the impossibility to perform the vertex analysis only, the cost function over the entire uncertain domain was calculated and then it will be shown for a more complete overview before the conventional costs vs. exibility index plots. It is worth remarking that the cost function associated to a given exibility index is the highest value over the domain covered by the corresponding deviation magnitude or, inversely, the index corresponds to the maximum deviation allowed for a given cost. Therefore, after minimizing the travel path for the system under perturbed conditions, the most conservative cost value should be accounted for.

Among the several randomly generated cases that were run, a common behaviour of the optimal solution and the related cost function could be detected.

In order to emphasize the results characteristic features, the most representative case will be used for the results discussion. It corresponds to the set of 10 cities whose TSP solution was already presented in Figure 1 with the addition of the starting depot vertex at the exact centre of the 2D space.

One run including the enumeration solution algorithm, the non-vertex deterministic exibility assessment and the stochastic one had an average computational time of about 6 minutes.

Distance uncertainty

The uncertainty concerning the distance between cities was studied rst. Several runs with randomly generated points and loads were performed in order to assess

# 1 0 H E A I C B J D F G 0 # 2 0 E I H C B J D F G A 0 # 3 0 E H C I A G F D J B 0 # 4 0 H C I E A G F D J B 0 # 5 0 H E I C B J D F G A 0
Tab. 2: Optimal paths -uncertain distances the critical aspects and select a case study worth to be presented. In general, the most relevant impact on the objective cost function was systematically observed for specic uncertain variables listed here below:

• the position of the closest city;

• the position of the farthest city;

• the position of the closest city with the heaviest load, i.e. highest l i /t 0i ratio;

• the position of the farthest city with the lightest load, i.e. lowest l i /t 0i ratio.

This phenomenon is due to the fact that they strongly aect the choice of the load to be delivered at the beginning or at the end of the travel by substantially modifying the optimal path.

An interesting set of points, as already anticipated, is the one whose TSP solution has been shown in Figure 1. Even in this case, the location of every point was perturbed and the cities who are farthest (B, A) and the closest (H, E) to the depot vertex resulted as the most sensitive ones as expected. In order to show the results in a 3D space, the two most relevant distances, namely 0B and 0H, have then been selected as uncertain variables. The expected deviation, represented by the dashed bars in Figure 1, is xed at 1.5 squares in the x vs y space since, for higher values, the travel cost function trend doesn't show any interesting behaviour but just keeps increasing linearly. Therefore, the exibility analysis for this case study was successfully carried out. Before

showing the exibility indices trends, the optimal TDP solutions and the value of the travel cost function over the uncertain, represented respectively in Figure 5a and Figure 5b, are worth to be discussed. Figure 5a shows that, despite the 362880 possible paths, only ve of them minimize the travel costs each one under specic perturbed conditions. It can be also noticed that the main dierence between them is the choice of the closest city to visit rst between H and E or the farthest one to visit between B/C and A/G on the top and bottom space regions respectively. Due to the overlapping of some paths, the list of the ve paths is provided in Table 2 for a better understanding.

On the other hand, the traveling cost function increases as the cities get farther, as expected. However, there are two more relevant aspects worth to be discussed in deepth. The rst one is the non-convexity of the cost function over the uncertain domain. This property is the reason why the exibility assessment performed by means of the simplied vertex analysis approach is not allowed for this specic problem. The second one is the position of the travel cost maximum for a given deviation magnitude view from the top; as for the previous Figures, the colormap represents the function value. In fact, as shown by the dotted line in Figure 5b, discontinuities in the trend of maximum costs can be observed. Each of them represents an optimal traveling sequence modication. In fact, while the optimal path can also change without showing such a trend discontinuity for a particular shape of the cost function, discontinuities, when present, are always related to an optimal path change.

This particular behaviour can be observed as well in the deterministic and stochastic exibility indices trends shown in Figures 6a and6b respectively.

The F SG index shows that, for the expected deviation, a 5.2 % additional travel costs under uncertain conditions are required. The index trend is almost linear with respect to the deviation magnitude and shows changes in its slope as discussed here above.

On the other hand, even if less remarkable, the SF index trend shows such discontinuities in its derivative as well both with the Beta and the Normal PDF. Dierently from the deterministic index, the stochastic one shows a nonzero value for the nominal operating conditions since the related travel costs already overestimates those obtained in other regions of the uncertain space as commented in others studies as well [START_REF] Di Pretoro | Flexibility analysis of a distillation column: Indexes comparison and economic assessment[END_REF]. The SF index with the Beta PDF starts from a lower value since the points with higher probabilities are displaced with respect to the nominal one but, for the same reason, it grows faster as the deviation magnitude increases. On the contrary, since the Gaussian PDF has the highest likelihood region centered in the nominal operating conditions, the related stochastic exibility index grows faster for low additional travel costs but is poorly improved for withstood deviations with higher magnitude. In conclusion, it can be stated that their behaviour reects that of the cumulative distribution functions expressed as a function of the deviation magnitude.

Load uncertainty

When it came to the delivery load uncertainty, the sensitivity analysis has been repeated on all the points in order to assess the most critical ones. Even in this case, they could be listed according to the four classes already presented in the previous section. For the specic values corresponding to the results presented here below, the most critical points were identied as the farthest point from the depot vertex, i.e. city B, and that with the highest l i /t 0i ratio, that is city E.

In this chapter we will show the results for a load deviation of both the points varying by the 90 % both in positive and negative directions since no signicant behaviour changes could be observed for higher value.

Figure 7a presents the optimal paths found during the exibility analysis.

Even in this case they are ve but, while four of them (2-5) are the same as those obtained before, the blue path, i.e. Path 1, has changed. In fact, for the particular set of loads, the point G has been replaced by the point I as a possible last city to visit before going back to the depot for very low values of l E . As in the previous section, due to the overlapping of some paths, the list of the ve paths is provided in Table 3 for a better understanding.

Even for this case, Figure 7b shows how the maximum costs trajectory varies with the deviation magnitude. In this case, the deviation refers to the loads to

# 1 0 H C B J D F G A E I 0 # 2 0 E I H C B J D F G A 0 # 3 0 E H C I A G F D J B 0 # 4 0 H C I E A G F D J B 0 # 5 0 H E I C B J D F G A 0
Tab. 3: Optimal paths -uncertain loads be delivered to the cities E and B, as it can be noticed on the two axes, and the value corresponds to the fraction with respect to its nominal value. Dierently from the previous section, although no direction changes can be detected, the travel cost function derivative shows a discontinuity (light green region) in the same direction. This behaviour can be noticed in the F SG trend for a 9 % deviation (cf Figure 8a).

In general, it can be observed that the deviation with respect to the delivery load has a higher impact on the traveling costs with respect to the city distance.

On the one hand this is due to the considerable uncertainty range taken into account while, on the other hand, the most relevant impact is due to the fact that a higher weight should be carried out during the entire travel that precedes the corresponding city.

As already commented, the Swaney and Grossman exibility index shows its usual almost linear trend with respect to the deviation magnitude with the exception of the slope discontinuity. The eect of this behaviour can be better observed in the SF graphics plotted in Figure 8b. This is explained by the fact that the stochastic exibility index accounts for data over the entire domain and not in correspondence of one point of the hyperrectangle perimeter only.

In fact, the Beta PDF shows an appreciable change in its convexity between 80

and 90 % and it approaches the unit value only when moving towards the top right (yellow) region of the uncertain domain. The same behaviour is much less relevant for the Gaussian PDF due to its symmetry and to the fact that the huge part of the perturbation likelihood lies in the proximity of small disturbances.

Conclusion

The initial purpose of this research work was to test the exibility assessment methodology and related indicators on an operational research problem and to point out analogies and dierences with respect to the conventional process application eld as well as to detect the eventual aspects to be adapted. In the light of the achieved result, it can be stated that the exibility indices available in literature t quite well this kind of problems with little adjustments required.

In general, this study shows that the adaption of the exibility assessment to the operational research domain allows to have more reliable expectations about the eective costs and the optimal operating solution variation in case external perturbations were likely to occur. Furthermore, this procedure makes possible the detection of the most critical supply chain parameters from a exibility perspective and their impact on the cost function so that further measures to be taken can be devoted to the system variables representing the main bottleneck. The key dierences to be noticed with respect to process equipment exibility analysis are the substantial discontinuity of the cost vs exibility trend slope and the more likely non-convexity of the feasible space. The reason of this behaviour was properly identied in the change of the optimal traveling sequencing when a certain disturbance magnitude is achieved.

Despite its simple formulation, the selected TDP case study properly shows the various critical aspects of supply chain management under uncertain operating conditions. In particular, the most critical parameters for the proposed example were identied into the position of the farthest and the closest points and in their loads respectively from a travel distance and cost function penalization perspective. The most important property of the employed methodology is its general validity. Indeed, normalized loads and traveling distances in a generic [-5 5]x[-5 5] y vs. x coordinates space have been used and a considerable number of randomly generated cases were run and always showed an analogous behaviour.

While for the deterministic exibility assessment, the deviation magnitude is the only piece of information required, for the stochastic one the impact of the uncertainty characterization, i.e. skewed vs symmetric PDF and the related parameter choice, was proved to be of critical importance both from a qualitative and quantitative perspectives. The higher information required by the SF is nevertheless well spent since, dierently from the deterministic indicator, its behaviour is not only the representation of the maximum expected costs trend but provides an overall perspective of the system performances.

For this reason, in case of non-convex objective functions, that is a common occurrence in MINLP optimization problems, the use of a stochastic indicator is particularly recommended.

Beside these relevant remarks, an additional important outcome of this study is that, given the operational research problem variables and the expected deviations, the optimal path under uncertain operating conditions always fall inside a limited set of alternatives despite the huge number of overall possible solutions. Moreover, these optimal traveling sequences dier from each other, and from the nominal solution, only by a small subsequence of cities according to the specic parameters that are subject to disturbances.

Due to the positive outcome of this research work, in future perspective, the proposed methodology is then worth to be extended to real industrial applications with more complex networks both from the supply and the demand perspective. For this purpose, dedicated software for analysis and processing of geospatial data from maps such as QGIS could be included to address problems of practical concern on existing systems. Furthermore, more eective algorithm from a computational point of view could be investigated in order to mitigate the computational time growth due to the higher complexity of the case study.

For instance, the same problem with more than two parameters at a time, such as the case of simultaneous distance and load deviations, could be easily addressed by applying the proposed procedure despite the impossibility to show the detailed results on a 3D uncertain space. Finally, the detailed economic assessment accounting both for the products and raw materials costs as well as logistic related ones could be included in order to have a complete overview about exible supply chain design.

Fig. 1 :

 1 Fig. 1: Traveling Salesman Problem solution

Fig. 2 :F

 2 Fig. 2: Deterministic exibility index F SG

Fig. 3 :

 3 Fig. 3: Stochastic exibility index SF

  the probability value related to each set of the uncertain parameters θ i lying inside the feasible region Ψ can be easily noticed for a general probability distribution function. The top view provided in the Figure shows the PDF values by means of a colormap from the lowest value (dark blue) to the highest (yellow).

  (a) Gaussian or Normal Distribution Function (b) Beta Distribution Function

Fig. 4 :

 4 Fig. 4: Bivariate Probability Distribution Functions

Fig. 5 :

 5 Fig. 5: Optimization results -uncertain distances

Fig. 6 :

 6 Fig. 6: Flexibility analysis results -uncertain distances

Fig. 7 :

 7 Fig. 7: Optimization results -uncertain loads

Fig. 8 :

 8 Fig. 8: Flexibility analysis results -uncertain loads