Flexible Assignment of Loading Bays for Efficient Vehicle Routing in Urban Last Mile Delivery

Urban freight deliveries are often subject to many access restrictions which creates the need to establish a system of loading bays and to split the last mile delivery into driving and walking parts. A new model based on hard and soft clustering approach is developed to solve the loading bay assignment problem for efficient vehicle routing and walking in last mile delivery. The flexibility of the model is provided by the soft clustering approach based on different membership degrees of customers to loading bays. Especially for instances with large numbers of loading bays, soft clustering seems to give better results, it leads to higher flexibility of city logistics systems, minimal driving distances, and adequately short walking paths, which contribute to the goal of reaching sustainable urban freight deliveries.


Introduction
In Europe, 72% of people live in cities, and this percentage is expected to rise to 80% by 2050 [1]. The importance of cities is also growing from a commercial point of view. Cities today already account for 85% of the EU's total gross domestic product (GDP) [2]. As a result, cities are experiencing a growing demand for transport with increasing congestion, noise, emissions, and other negative effects of transport on the environment and city dwellers [3,4]. These problems are most prevalent in historical cities and city centers, which are generally not designed to cope with such high traffic volumes. These urban areas have major accessibility problems (narrow roads) and a lack of logistics infrastructure [5]. Cities are further hampered by the modern trends of just-in-time delivery and e-commerce, which increase the number of deliveries and consequently lead to fragmentation of urban freight transport [6,7]. This leads to an increase in the number of journeys and low utilization of delivery vehicles [8]. Studies show that lorries account for 15% of urban traffic and more than 20% of CO 2 emissions and congestion [9].
Limited space in urban centers does not allow for further expansion of transport and logistics infrastructure. Cities are therefore often restricting access of delivery vehicles to urban centers and see this approach as the only option to ensure sustainable living conditions [10]. The restrictions usually relate to the time of the day and the emission class or size of the delivery vehicles. Strict city centers are often defined as pedestrian zones where access conditions are further restricted. Despite the best intentions of city authorities, such measures worsen the accessibility of locations in city centers, negatively influencing attractiveness for commercial services [11].
found in the literature [18,28,29] and to the best of our knowledge, there are no examples of 2E-LRP dealing with loading bays and walking parts of urban last mile delivery.
Scholars working on urban freight transport issues try to solve these complex problems mainly by using specially developed and tailored simulation-based models [30,31]. In the literature we have found only four models that specifically address the cases of loading bays. The first model was developed by Comi et al. [32] who used the discrete event simulation approach to investigate the possibilities of loading bay reservation for last mile delivery in the limited traffic zone of Rome. The second is provided by Zhang and Thompson [33], who developed an agent-based model to study the importance of providing loading bays occupancy information on driving and walking costs. The third model was developed by Mc Load et al. [34], who used a metaheuristic approach to quantify the environmental and financial benefits of using porter and cycle couriers for specific last mile parcel delivery in London. The last one we found was by Nguyen et al. [35], who developed a mathematical model to optimize total delivery time based on dual-mode route planning. The second mode relates to walking time of deliverers and considers them as an important element of the vehicle routing procedure. The main feature of the previously mentioned models is that they all, except Zhang and Thompson [33], deal with a specific real-life problem that results in a very specific solution. This approach is of course very efficient for this specific example, but it does not allow to understand basic (general) principles that could be used for planning loading bay facilities under different circumstances. The other important issue relates to the flexibility of such a system and the dynamic assignment of loading bays if the best loading bay is occupied.
This article attempts to fill this gap by developing a model based on flexible assignment of loading bays for the vehicle routing in urban last mile delivery. The general research aim of the article is not to find the optimal solution for both echelons (routing and walking), but to investigate the efficiency of a flexible assignment of loading bays for vehicle routing and consequently for the walking part of last mile delivery. The solution to this problem is very important, as solutions would enable city authorities to better manage urban freight deliveries and the results could be applicable, to a certain extent, to similar problems in other research areas. For example, the problem of assigning incoming vehicles to loading bays is mathematically similar to the berth allocation problem of ships in seaports [36][37][38], the problem of locating urban public e-charging stations [39], the location of various public services [40], and different types of collection points [41].

Materials and Methods
In Letnik et al. [42] a multi-parametric model of the idealized urban area was introduced and applied to different instances of rectangular urban grid structured zones to determine statistically the most relevant number and location of loading bays. The main criterion for the evaluation of the results was the predefined maximum allowed walking distance from the loading bay to the customer. In this paper the existing model is upgraded with vehicle routing algorithm to address also this part of the delivery. The main idea is to investigate the flexible loading bay assignment procedure when the delivery is made to a specific customer or a group of customers in the city center. Many fundamentally different cases can occur where several loading bays are available for delivery to a particular customer and where more customers need to be served from the same loading bay.
To solve this issue, we developed a methodology and examined two different assignment methods to understand which approach is better from the driving and walking distance point of view. In the first approach, only the closest loading bay for each customer in the city center is selected, whereas in the second one, up to three potential loading bays (using the soft clustering approach) can be selected and used. In the continuation of this paper, the first approach is referred to as "hard" and the second as "soft," which implies to the hard and soft clustering procedure. The model is based on the concept of the idealized urban area, which allows us to simulate and understand basic principles regardless of the characteristics of a particular real-life case. The above-described problem was broken down into two separate sub-problems, although we are aware that both sub-problems (vehicle routing and on foot delivery) are interconnected. Consequently, the separate solutions of the two sub-problems would therefore at best be an approximation of the optimal solution of the overall problem. Due to approximation, on the other hand, such a simplified approach enables a much more efficient solution to the described problem and of similar real-life problems. Moreover, all assumptions considered in the presented model are based on the data from the cited literature and/or on our own experience, and they follow scientific principles that allow us to easily understand, explain, and generalize the problems under consideration.

Model Characteristics and Functions
The main elements of the model are the road transport network, access gates, customers' locations, and location of loading bays (see Figure 1). Idealized urban area is presented graphically as an orthogonal network of roads on which customers, loading bays, and access gates are located. Customers are placed on the road network randomly, following the uniform density. For this purpose, Latin hypercube sampling (LHS) algorithm is used [43]. The positions of loading bays are determined by the fuzzy C-means (FCM) procedure [42], access gates are located randomly on each side of the urban area perimeter also using LHS algorithm.
Sustainability 2020, 12, x FOR PEER REVIEW 4 of 19 principles that allow us to easily understand, explain, and generalize the problems under consideration.

Model Characteristics and Functions
The main elements of the model are the road transport network, access gates, customers' locations, and location of loading bays (see Figure 1). Idealized urban area is presented graphically as an orthogonal network of roads on which customers, loading bays, and access gates are located. Customers are placed on the road network randomly, following the uniform density. For this purpose, Latin hypercube sampling (LHS) algorithm is used [43]. The positions of loading bays are determined by the fuzzy C-means (FCM) procedure [42], access gates are located randomly on each side of the urban area perimeter also using LHS algorithm. Note that the model allows to create different instances of orthogonal urban area with different density of network, customers, and loading bays. Example of an idealized city plan with random location of customers (green circles), optimal position of loading bays (red dots), and random location of access gates (black squares) are presented in the Figure 2 below. As can be seen, we assume one can enter and exit the urban area only at specific locations (access gates). This is a reasonable assumption because this is often the case in urban areas with implemented regulated urban access restriction. Note that the model allows to create different instances of orthogonal urban area with different density of network, customers, and loading bays. Example of an idealized city plan with random location of customers (green circles), optimal position of loading bays (red dots), and random location of access gates (black squares) are presented in the Figure 2 below. As can be seen, we assume one can enter and exit the urban area only at specific locations (access gates). This is a reasonable assumption because this is often the case in urban areas with implemented regulated urban access restriction. density of network, customers, and loading bays. Example of an idealized city plan with random location of customers (green circles), optimal position of loading bays (red dots), and random location of access gates (black squares) are presented in the Figure 2 below. As can be seen, we assume one can enter and exit the urban area only at specific locations (access gates). This is a reasonable assumption because this is often the case in urban areas with implemented regulated urban access restriction.  The model contains two key optimization modules: FCM module, which generates a set of customers locations and provides the input for multiple facility location problem (MFLP), in our case multiple locations of loading bays. Vehicle routing module, which takes as a basis the results of FCM, selects optimal loading bays for specific customers while providing the solution for optimal vehicle routing. FCM is based on the principles of fuzzy logic, while the mathematical origin of vehicle routing is based on integer linear programming. Both procedures were implemented in MATLAB and are briefly presented below.

Clustering Procedure for Solving Facility Location Problem
Multiple facility location problem can be treated as a classical clustering problem, where customers and facility locations correspond to data points and centers. The clustering criterion is to minimize the sum of distances from each data point to the corresponding cluster center. Potential location of loading bays in the model is therefore determined by using fuzzy C-means (FCM) clustering procedure. FCM was chosen to get some flexibility to the model with the possibility to simulate different membership degrees of customers to cluster centers, in our case loading bays [42].
FCM aims at minimizing the following function [44]: where: -D = number of customers -N = number of clusters - x i = i-th data point (each data point corresponds with location of a particular customer) -c j = center of cluster j µ ij = degree of membership of customer i to cluster j m = parameter of fuzziness FCM algorithm in implemented in MATLAB and performs the following steps: 1.
Randomly initialize the cluster membership values, µ ij .
Update µ ij according to the following: Calculate the objective function, J m . 3.
Repeat steps 2-4 until J m improves by less than a specified minimum threshold ε or until after a specified maximum number of iterations.

Vehicle Routing and Loading Bay Assignment Procedure
Urban last mile delivery mostly involves three shipments while the number of packages per delivery varies significantly depending on the type of cargo and the type of customers [45,46]. Based on this assumption we developed a model which simulates three shipments per delivery which is called a triple delivery. The essence of the simulations is to find the optimal vehicle routing for triple delivery according to the randomly selected three customers C 1 , C 2 , C 3 , the corresponding loading bays L 1 , L 2 , or L 3 , and access gates A 1 , A 2 , A 3 , or A 4 . Schematic representation of a graph is presented in the Figure 3. Urban last mile delivery mostly involves three shipments while the number of packages per delivery varies significantly depending on the type of cargo and the type of customers [45,46]. Based on this assumption we developed a model which simulates three shipments per delivery which is called a triple delivery. The essence of the simulations is to find the optimal vehicle routing for triple delivery according to the randomly selected three customers , , , the corresponding loading bays , , or , and access gates , , , or . Schematic representation of a graph is presented in the Figure 3. Blue edges in the graph ( Figure 3) represent connections between access gates and loading bays ( , ), orange edges represent connections among loading bays ( , ), and green edges represent connections between loading bays and customers ( , ). The customers are not connected because we assume all deliveries are made individually even if one loading bay is selected for several customers.
In order to solve the problem we first consider the hard optimization approach at three loading bays , , , and then also the soft optimization approach on alternative loading bays , , , . Loading bays are assigned to customers based on FCM clustering procedure considering different membership degrees. Depending on the position of the selected customers and the associated loading bays, the solution of vehicle routing can be, in case of using soft optimization approach, drastically improved, especially if = for any , , , ∈ {1,2,3}. It is very important to clarify that the aim of supplying the city customers through the limited number of access gates is to avoid "unnecessary" traffic within the urban area. From that point of view, we assumed that the distances between individual access gates are irrelevant. In practice, this Blue edges in the graph ( Figure 3) represent connections between access gates and loading bays (x i,j ), orange edges represent connections among loading bays (x j,j ), and green edges represent connections between loading bays and customers (x j,c ). The customers are not connected because we assume all deliveries are made individually even if one loading bay is selected for several customers.
In order to solve the problem we first consider the hard optimization approach at three loading bays L 1 , L 2 , L 3 , and then also the soft optimization approach on alternative loading bays L 11 , L 12 , L 13 , L 21 , L 22 , L 23 , L 31 , L 32 , L 33 . Loading bays are assigned to customers based on FCM clustering procedure considering different membership degrees. Depending on the position of the selected customers and Sustainability 2020, 12, 7500 7 of 19 the associated loading bays, the solution of vehicle routing can be, in case of using soft optimization approach, drastically improved, especially if L ij = L kl for any i, j, k, l ∈ {1, 2, 3}.
It is very important to clarify that the aim of supplying the city customers through the limited number of access gates is to avoid "unnecessary" traffic within the urban area. From that point of view, we assumed that the distances between individual access gates are irrelevant. In practice, this means that the costs, and therefore the corresponding negative effects of moving the vehicle from one access gate to another, can be neglected. This is the main assumption and limitation of the described model (of course, the model also allows modelling on the assumption that the distances between the gates are equal to the distances around the perimeter of the square or rectangle on which we are modelling).
The aim of the model is therefore to solve the problem of vehicle routing together with choosing a combination of loading bays and access gates (related to blue and orange lines in Figure 3) and to minimize the sum of all driven kilometers within the city center. The minimum of the objective function (i.e., kilometers) is sought by means of combinatoric optimization via breaking down the problem into all possible outcomes and apply integer linear programming to each outcome in order to find the overall minimum of driven kilometers.
In order to handle the routing integrated with location-allocation problem, there are 54 types of different vehicle routing outcomes. For each outcome, a unique linear program provides the partial optimal solution. The overall optimal solution is the minimum of all partial optimal solutions. Below, we present the necessary inputs to formulate these integer linear programs and explicitly define one of the outcomes and the corresponding integer linear program.
First we need to calculate distances between access gates and loading bays (cA i L j ) and between loading bays (cL i L j ; i j). Distances are calculated using Manhattan metrics in the following way: and In the above equations abs stands for the absolute value and (X, Y) is used for the physical point where the indices refer to the particular physical point (e.g., Y L j is the Y coordinate of the point L j ).
All possible combinations (calculated distances) are written as a cost vector: Furthermore, the point vector is denoted by: Next, define rows in the matrix: which presents the constraints of the integer linear program. For example, to allow that the gate A 1 is connected to one of the loading bays L 1 , L 2 or L 3 we define the following line vector: In a similar way we define {mA 2 }, {mA 3 }, and {mA 4 } = {0 0 0 1 0 0 0 1 0 0 0 1 0 0 0} for allowing that the other three gates are connected to one of the loading bays.
To connect L 1 with some gate A i and with some another loading bay; L 2 or L 3 we define the following line vector: and similarly we define {mL 2 } and {mL 3 } = {0 0 0 0 0 0 0 0 1 1 1 1 0 1 1}. Finally, the right-hand-side vector, {rhs}, of the matrix-form constraint equation must be defined. For example, if access gates A 1 and A 3 are used and the delivery is performed through all three loading bays, we have Using these notations, the integer linear program for solving optimization of the presented problem, writes as follows: For example, for the case of using assess gates A 1 and A 2 and supplying customers each from its own loading bay; L 1 , L 2 , and L 3 , the optimal solution can be found with the help of linear program aiming to minimize the following cost function: subject to x 1,5 + x 1,6 + x 1,7 = 1, where for all i, j ∈ {1, 7} For clarification, the constraint x 1,5 + x 1,6 + x 1,7 = 1 means that the solution is starting from or ending in the gate A 1 . The constraint x 2,5 + x 2,6 + x 2,7 = 0 means that the gate A 2 is not used in the solution. The constraint x 1,5 + x 2,5 + x 3,5 + x 4,5 + x 5,6 + x 5,7 = 2 means that the solution enters and leaves loading bay L 1 .
In a similar way we can treat all other feasible solutions for triple delivery using three, two, or one loading bay.

Potential Solutions/Outcomes
To list all possible outcomes, we use abbreviations which uniquely determine any combination of access gates and loading bays. The first part of the notation defines the number of gates (G#), the second part the number of loading bays (L#), and the final part the type of the solution, where: P stands for a path and C for a cycle. The notation PAiA j stands for a path from Ai to A j, CnAi stands for a cycle on n points, n ∈ {2, 3, 4}, containing gate Ai. And finally, EAi denotes the gate Ai is empty, which becomes relevant only if the customers are each supplied through its own entrance gate.
We are first listing four significantly different outcomes for the case of supplying customers through three loading bays. The first example is to check the solution, where we enter the city through three different access gates A1, A2, and A3 and supply customers each from different loading bays L 1 , L 2 , and L 3 . In this case the entire route can be formally written as: G3-L3-EA4. Another example is to check the outcome, where we enter the city through the access gate A 1 and exit through gate A 2 and we supply customers each from different loading bays; L 1 , L 2 , and L 3 . The entire route can be then formally written as: G2-L3-PA1A2. Next outcome is describes the situation when we enter the city Sustainability 2020, 12, 7500 9 of 19 through the access gate A 1 and exit through gate A 3 to serve two customers each from different loading bay and then to supply the third customer entering and exiting the city through gate A 3 . The entire route can be in this case formally written as: G3-L3-PA1A2C2A3. The final example describes the outcome, where we enter the city through only one access gate, A 1 , and we supply customers each from its own loading bay; L 1 , L 2 , and L 3 . In this final case the entire route can be formally written as: G1-L3-C4A1. To make clearer, the above outcomes are graphically presented in Figure 4.
We are first listing four significantly different outcomes for the case of supplying customers through three loading bays. The first example is to check the solution, where we enter the city through three different access gates A1, A2, and A3 and supply customers each from different loading bays , , and . In this case the entire route can be formally written as: G3-L3-EA4. Another example is to check the outcome, where we enter the city through the access gate and exit through gate and we supply customers each from different loading bays; , , and . The entire route can be then formally written as: G2-L3-PA1A2. Next outcome is describes the situation when we enter the city through the access gate and exit through gate to serve two customers each from different loading bay and then to supply the third customer entering and exiting the city through gate . The entire route can be in this case formally written as: G3-L3-PA1A2C2A3. The final example describes the outcome, where we enter the city through only one access gate, , and we supply customers each from its own loading bay; , , and . In this final case the entire route can be formally written as: G1-L3-C4A1. To make clearer, the above outcomes are graphically presented in Figure 4. In the case of using three loading bays we must in total execute 38 different integer linear programs with different right-hand-side vectors {rhs}.
Next option is to supply customers through only two loading bays. This means that two out of three customers are served from the same loading bay. We present three examples of significantly different outcomes for this option. The first example is to check the solution, where we enter the city through two different access gates, A 1 and A 2 . From the gate A 1 we visit the first loading bay L 1 to supply the first customer and from gate A 2 we visit second loading bay L 2 to supply the other two customers. In both cases we exit the city through the same gate as we entered. This can be formally written as: G2-L2-C2A1C2A2. Next example is to check the solution of entering the city through access gate A 1 and exiting through gate A 2 . On the path we visit the first loading bay L 1 to supply two customers and then continue to second loading bay L 2 to supply the remaining customer. This can be formally written as: G2-L2-PA1A2. Last example is to check the solution, where we enter the city through only one access gate A 1 and we supply customers from two loading bays L 1 and L 2 . From the first loading bay we supply two customers and from the second the remaining one. The entire route can be in this case formally written as: G1-L2-C3A1. For easier understanding, the before-mentioned examples are also graphically and schematically presented in Figure 5.  In the case of using three loading bays we must in total execute 38 different integer linear programs with different right-hand-side vectors { ℎ }.
Next option is to supply customers through only two loading bays. This means that two out of three customers are served from the same loading bay. We present three examples of significantly different outcomes for this option. The first example is to check the solution, where we enter the city through two different access gates, and . From the gate we visit the first loading bay to supply the first customer and from gate we visit second loading bay to supply the other two customers. In both cases we exit the city through the same gate as we entered. This can be formally written as: G2-L2-C2A1C2A2. Next example is to check the solution of entering the city through access gate and exiting through gate . On the path we visit the first loading bay to supply two customers and then continue to second loading bay to supply the remaining customer. This can be formally written as: G2-L2-PA1A2. Last example is to check the solution, where we enter the city through only one access gate and we supply customers from two loading bays and . From the first loading bay we supply two customers and from the second the remaining one. The entire route can be in this case formally written as: G1-L2-C3A1. For easier understanding, the beforementioned examples are also graphically and schematically presented in Figure 5. In case when using two loading bays we must execute 16 different integer linear programs with different right-hand-side vector { ℎ }.
The last option is to supply three customers simultaneously through only one loading bay. This means that customers are located at a relatively short distance from each other, so the best solution is to stop the truck at one loading bay and supply all customers from there. In this case, only the minimum distance from the selected loading bay to the nearest access gate needs to be found. Therefore, we do not need to run a linear program, we just need to search for the shortest path. In case when using two loading bays we must execute 16 different integer linear programs with different right-hand-side vector {rhs}. The last option is to supply three customers simultaneously through only one loading bay. This means that customers are located at a relatively short distance from each other, so the best solution is to stop the truck at one loading bay and supply all customers from there. In this case, only the minimum distance from the selected loading bay to the nearest access gate needs to be found. Therefore, we do not need to run a linear program, we just need to search for the shortest path.

Simulation Results
The aim of the simulations is to perform vehicle routing optimization for specific instances of idealized urban areas in order to compare the results of driving and walking distances for different numbers of loading bays using soft and hard optimization approach and to get, if they exist, general relations/rules for vehicle routing phenomenon. Simulations are carried out on a square grid with density of ρ n = 11 × 11 roads, customer density of ρ c = 1000 customers/km 2 , and for the following number of loading bays, #L = 4, 9, 16, 25, 49, and 81. These parameters are carefully selected to represent a typical urban area of European cities as well as possible. For more details see Letnik et al. [47].
The following Figure 6 show some basic examples of idealized urban areas with randomly generated customer locations, corresponding locations of loading bays (generated by FCM procedure), and randomly located access gates. These "instances" are examples of different situations that can occur in real-life. Analyzing large numbers of such instances is the cornerstone of our statistical approach, which first serves to locate the loading bays and then to understand the influence of a flexible assignment of loading bays for vehicle routing and consequently walking distances. Note that these instances form the basis for the simulations shown below. statistical approach, which first serves to locate the loading bays and then to understand the influence of a flexible assignment of loading bays for vehicle routing and consequently walking distances. Note that these instances form the basis for the simulations shown below. Vehicle routings and the corresponding walking parts are calculated separately for each instance presented in Figure 6. The algorithm randomly selects three customers for a given instance and then executes all together 54 different types of integer linear programs to find an optimal solution for the driving part to serve these three customers. When the solution of the integer linear program is obtained, the algorithm selects the next triple. The simulations are performed for the density of =1000 customers/ , where 333 simulations of a triple delivery are performed for each instance. Note that in the case of hard (soft) optimization approach, the algorithm considers up to three (nine) different potential loading bays for each simulation of triple delivery. This yields that in the statistical analysis we considered 333 • 54 • 2 simulations for each instance. It is important to mention that optimum is sought from the driving distance point of view, whereby walking part is kept reasonably low by assigning only those loading bays with reasonably high membership values.
Example of expected result of vehicle routing in case of hard and soft clustering approach is presented in Figure 7. Vehicle routings and the corresponding walking parts are calculated separately for each instance presented in Figure 6. The algorithm randomly selects three customers for a given instance and then executes all together 54 different types of integer linear programs to find an optimal solution for the driving part to serve these three customers. When the solution of the integer linear program is obtained, the algorithm selects the next triple. The simulations are performed for the density of ρ c = 1000 customers/km 2 , where 333 simulations of a triple delivery are performed for each instance. Note that in the case of hard (soft) optimization approach, the algorithm considers up to three (nine) different potential loading bays for each simulation of triple delivery. This yields that in the statistical analysis we considered 333·54·2 simulations for each instance. It is important to mention that optimum is sought from the driving distance point of view, whereby walking part is kept reasonably low by assigning only those loading bays with reasonably high membership values.
Example of expected result of vehicle routing in case of hard and soft clustering approach is presented in Figure 7.
Vehicle routings and the corresponding walking parts are calculated separately for each instance presented in Figure 6. The algorithm randomly selects three customers for a given instance and then executes all together 54 different types of integer linear programs to find an optimal solution for the driving part to serve these three customers. When the solution of the integer linear program is obtained, the algorithm selects the next triple. The simulations are performed for the density of =1000 customers/ , where 333 simulations of a triple delivery are performed for each instance. Note that in the case of hard (soft) optimization approach, the algorithm considers up to three (nine) different potential loading bays for each simulation of triple delivery. This yields that in the statistical analysis we considered 333 • 54 • 2 simulations for each instance. It is important to mention that optimum is sought from the driving distance point of view, whereby walking part is kept reasonably low by assigning only those loading bays with reasonably high membership values.
Example of expected result of vehicle routing in case of hard and soft clustering approach is presented in Figure 7. As can be seen in Figure 7, the hard clustering approach takes into consideration only loading bays with the highest membership degree (L 11 , L 21 , L 31 ), while the soft clustering approach considers also loading bays with the second and third highest membership degrees (L 12 , L 13 , L 22 , L 23 , L 32 , L 33 ). Consequently, shortest driving distances are expected in case of soft clustering approach, as can be seen from the example presented above. Note that on the other hand this results in slightly longer walking distances.
The comparison between the driving and walking distances in the case of a hard and a soft optimization approach is in the following analyzed and presented by cumulative distribution function (CDF) curves. We chose this method because it allows to read simultaneously the minimum, mean, and maximum distances on the same graph, and therefore we can easily compare the differences between hard and soft optimization approach.

Results-Driving Distances
The understanding of driving distances is important to assess the impact of delivery processes on the environmental parameters of the city. The aim of optimizing delivery is namely to reduce, as much as possible, the kilometers driven within the urban area. In the analysis we are trying to understand the influence of the number of loading bays on the driving distances and their distribution in the overall structure of the delivery. In the following Figures 8-15 and Tables 1-4 "HD" stands for hard driving and "SD" stands for soft driving solution. As can be seen in Figure 7, the hard clustering approach takes into consideration only loading bays with the highest membership degree ( , , ), while the soft clustering approach considers also loading bays with the second and third highest membership degrees ( , , , , , ). Consequently, shortest driving distances are expected in case of soft clustering approach, as can be seen from the example presented above. Note that on the other hand this results in slightly longer walking distances.
The comparison between the driving and walking distances in the case of a hard and a soft optimization approach is in the following analyzed and presented by cumulative distribution function (CDF) curves. We chose this method because it allows to read simultaneously the minimum, mean, and maximum distances on the same graph, and therefore we can easily compare the differences between hard and soft optimization approach.

Results-Driving Distances
The understanding of driving distances is important to assess the impact of delivery processes on the environmental parameters of the city. The aim of optimizing delivery is namely to reduce, as much as possible, the kilometers driven within the urban area. In the analysis we are trying to understand the influence of the number of loading bays on the driving distances and their distribution in the overall structure of the delivery. In the following Figures 8-15 and Tables 1-4 "HD" stands for hard driving and "SD" stands for soft driving solution.   As can be seen, driving distances are in all cases shorter when using soft optimization approach. It can also be observed that the difference in driving distances between hard and soft optimization approach decreases with the increase in the number of loading bays.
The differences are shown in more detail in the following Table 1.

Results-Walking Distances
Walking distances represent the second part of urban last mile delivery. In comparison to driving distances, walking does not cause any negative effects on the environment but still needs to be as short as possible to avoid longer stay of delivery vehicles in the urban area. Also in this case we need to understand the differences of hard and soft optimization approach. In the following figures As can be seen, driving distances are in all cases shorter when using soft optimization approach. It can also be observed that the difference in driving distances between hard and soft optimization approach decreases with the increase in the number of loading bays.
The differences are shown in more detail in the following Table 1.

Results-Walking Distances
Walking distances represent the second part of urban last mile delivery. In comparison to driving distances, walking does not cause any negative effects on the environment but still needs to be as short Sustainability 2020, 12, 7500 13 of 19 as possible to avoid longer stay of delivery vehicles in the urban area. Also in this case we need to understand the differences of hard and soft optimization approach. In the following figures and tables "HW" stands for hard walking and "SW" stands for soft walking solution.
Sustainability 2020, 12, x FOR PEER REVIEW 13 of 19 Figure 11. CDF for walking distances, hard and soft approach, # = 4 and # = 9  As can be seen, walking distances are longer when using soft optimization approach which is exactly the opposite as for the case of driving distances presented before. It can also be observed that the difference between hard and soft approach decreases with the increase in the number of loading bays which is similar as by the case of driving distances.
The differences are shown in more detail in the following Table 2.  As can be seen, walking distances are longer when using soft optimization approach which is exactly the opposite as for the case of driving distances presented before. It can also be observed that the difference between hard and soft approach decreases with the increase in the number of loading bays which is similar as by the case of driving distances.
The differences are shown in more detail in the following Table 2.  As can be seen, walking distances are longer when using soft optimization approach which is exactly the opposite as for the case of driving distances presented before. It can also be observed that the difference between hard and soft approach decreases with the increase in the number of loading bays which is similar as by the case of driving distances.
The differences are shown in more detail in the following Table 2.  As can be seen, walking distances are longer when using soft optimization approach which is exactly the opposite as for the case of driving distances presented before. It can also be observed that the difference between hard and soft approach decreases with the increase in the number of loading bays which is similar as by the case of driving distances.
The differences are shown in more detail in the following Table 2. As can be seen from Table 2, walking distances are, in the case of using soft optimization approach, almost exactly double of those reached by hard optimization approach. In addition, this holds true for different number of loading bays.

Comparison Between Driving and Walking Distances
Our final aim is to understand the relationship and compare the differences between the driving and walking distances by different number of loading bays. First, we show a comparison for the hard optimization approach. As can be seen from Table 2, walking distances are, in the case of using soft optimization approach, almost exactly double of those reached by hard optimization approach. In addition, this holds true for different number of loading bays.

Comparison Between Driving and Walking Distances
Our final aim is to understand the relationship and compare the differences between the driving and walking distances by different number of loading bays. First, we show a comparison for the hard optimization approach. From the chart above we can see:

•
The total length of average driving and walking distances decreases with the increase in the number of loading bays. • The driving distance covered in the hard optimization approach does not show significant deviations with regard to the increase in the number of loading bays.
• Increasing the number of loading bays on the other hand shortens the length of walking distances.
Distances and differences are shown in more detail in the following Table 3.  From the chart above we can see: • The total length of average driving and walking distances decreases with the increase in the number of loading bays.

•
The driving distance covered in the hard optimization approach does not show significant deviations with regard to the increase in the number of loading bays.

•
Increasing the number of loading bays on the other hand shortens the length of walking distances.
Distances and differences are shown in more detail in the following Table 3.  We can also show a comparison for the soft optimization approach. From the chart above we can see:

•
The total distance (traveling and walking) decreases with the increase in the number of loading bays (similar to what we noticed with the hard optimization approach).

•
The ratio between the length of the driving and walking distances is significantly different in this case; with the increase in the number of loading bays, there is a pronounced trend of increasing the length of driving distances and decreasing the length of walking (except for # = 16, the length of driving distance is slightly longer than for # = 49).
Distances and differences are shown in more detail in the following Table 4.  From the chart above we can see:

•
The total distance (traveling and walking) decreases with the increase in the number of loading bays (similar to what we noticed with the hard optimization approach).

•
The ratio between the length of the driving and walking distances is significantly different in this case; with the increase in the number of loading bays, there is a pronounced trend of increasing the length of driving distances and decreasing the length of walking (except for #L = 16, the length of driving distance is slightly longer than for #L = 49).
Distances and differences are shown in more detail in the following Table 4.

Discussion
Urban freight deliveries are often subject to many access restrictions which creates the need to establish a system of loading bays and to split the last mile delivery into driving and walking parts. An open question arises-where and how many loading bays should we provide to perform this transshipment and execute last mile delivery in the most efficient way.
The problem is even more complex because loading bays are sometimes occupied and, in these cases, we need to reroute delivery vehicles and search for the next alternative loading bay. To test different optimization approaches we introduced fuzzy clustering to make the system flexible enough for accommodating also this kind of problem. To stress once again, the approach is based on the following important assumptions:

•
The modelled city corresponds to idealized orthogonal urban area with different densities of loading bays; • Locations of customers and access gates are generated randomly, following the uniform density (simulated by LHS procedure); • Vehicles can enter and exit the urban area only at specific locations (access gates) and the distances between them (outside the studied area) are considered irrelevant; • Simulations are restricted on up to three shipments per delivery (triple delivery); • All deliveries from loading bays to customers are made individually even if one loading bay is selected for several customers.
The aim of the simulations was to understand the influence of different numbers of loading bays on driving and walking distances when using hard and soft (fuzzy) optimization approach. An important finding of the simulation analysis is related to the case of hard optimization approach. We show there is no statistically significant influence of the number of loading bays on the driving distances. In practical terms this means that city authorities are not able to decrease total driving distances with increasing the number of loading bays. The crucial decision parameter for determining the number of loading bays is therefore only the acceptable walking distance, which decreases with increasing numbers of loading bays. This is a very important finding of the model, which we would not encounter if we would run the simulation for specific characteristics of a particular urban area. This result can be seen as a fundamental law that is relevant for all rectangular urban areas.
On the other hand, soft optimization approach shows its efficiency for decreasing driving distances, in particular when the number of loading bays is small. Whereby, in this case, the walking part of the delivery is much higher than in the case of hard clustering optimization approach. The situation on the other hand changes with increasing numbers of loading bays where the walking distances for soft optimization approach becomes much lower and acceptable also from the practical point of view. Lower walking distances are also influencing the time of loading bay occupancy and are therefore additionally important, especially in cases where loading bays have some capacity restrictions. This result shows the importance of understanding the basic relation between walking and driving distances when using different optimization approaches.
We can conclude that a soft optimization approach proves to be an acceptable alternative for decreasing driving distances especially for a larger number of loading bays. We would therefore advise city authorities to plan reasonably large numbers of loading bays in urban areas and to use soft optimization approach when assigning loading bays to the routing of the urban last mile deliveries. This would enable them to reduce the number of driven kilometers and create greater flexibility of the city logistics system. Such a system would also have beneficial results for logistics operators, users, and city dwellers, which are constantly looking for a better, more efficient and sustainable city logistics system.
Before finally confirming this results on some real life cases, further research could also investigate the impact of different densities of the road network (including variants with an uneven distribution of roads that are closer to some specific real-life situations), different densities of customers, and different numbers and location of access gates. Nevertheless, on the results presented in this article and the research experience gained in determining the optimal number and locations of loading bays [42], we assume that the results will not be significantly different from those shown. Another proposal is to run the optimization also for the walking part of delivery, which would additionally save some walking distances.