Genetic Algorithm-Based Optimization for the Fuzzy Capacitated Location-Routing Problem with Simultaneous Pickup and Delivery

A significant portion of a company's budget is typically allocated to logistics costs. These expenses can be minimized through the meticulous design of the supply chain, which serves as a crucial distribution network. Ultimately, this network comprises numerous small product flows directed toward customers or final vendors, ensuring efficient delivery [1]. The design of the network leads to two hard combinatorial optimization problems: determining the location of warehouses and determining vehicle routes for giving service to customers of the warehouses. These two kinds of decisions used to be taken separately for a long time, but continued progress in optimization techniques makes it possible to integrate the approaches as LRP [2]. In fact, LRP aims at finding suitable locations, numbers of facilities, and distribution routes for vehicles.

The evolution trend of the logistics model of optimization can be summarized as follows. At the beginning, LRP classic models were examined separately. Afterwards, there was a need for integrating the two issues. As a result, a location-routing model was introduced. The problem concerned the cases where only customers received their orders, and they were not allowed to send goods and materials to warehouses and distribution centers. However, in practice, there are many cases where customers, in addition to receiving the goods, also send them, e.g., glass beverage distribution system, distribution of dairy products every day, etc., so the applied model of LRPSPD was discussed. However, the most challenging and problematic problem in all areas dealing with the real world is that in many cases, there is no certainty, and the demands of customers are often unknown. This means that the information relating to demand and pickup in each of the customers often is not precise enough. For example, based on experience, it can be concluded that customer demands (about 50 units) are between 20 and 60 units. For this reason, most of the used data are not precise enough for the probability distribution of customer demand. On the other hand, based on the judgement of experts, you can easily meet the demands of customers. Therefore, while the use of probability theory is difficult and costly, fuzzy logic can be useful in resolving this problem [3]. Accordingly, we can conclude that an LRPSPD model with a fuzzy approach to what extent could be close to reality and applied.

The concept of supply chain was first introduced in the late 1980s and extended later in the 1990s. This concept can be considered one of the most important trends in management and operation in the last two decades. One important reason for this claim is the high share of logistics costs in the cost of a product (about 30%). The supply chain includes all measures that directly or indirectly cause the needs of customers to be met. It does not include suppliers and builders but also includes transportation, facilities, retailers, and the customers themselves (Figure 1).

The location-allocation problem is one of the important problems as well as a location measure for modelling in logistics systems. The aim of solving this problem is the simultaneous determination of the optimal number and location of facilities and the allocation of customers to facilities, assuming delivery of customer demand by a full cargo vehicle. Cooper introduced and solved the location-allocation problem in 1963 [4].

The vehicle routing problem (VRP), which was first introduced by Dantzig and Ramser [5], is one of the most important optimization problems from research and operational perspectives. In this problem, we look to determine a set of routes, each of which is used by a vehicle, so that we meet the customer needs, satisfy all operational constraints, and minimise the total cost. In this problem, vehicles departing from a distribution centers return to the same distribution center.

Components of LRP include the location problem, the allocation problem, and the routing problem.

LRP has been one of the most important problems of location among the logistical systems analysts in the past three decades. These problems provide analysts and executives with a strong perspective and prevent local optimization of the LRP because of considering the interaction between facility location and routing of transportation. The main difference between LRP and the location-allocation problem is that after determining the location of the facility in the location-allocation problem, travels occur directly from the facilitator to the supplier/customer, but in LRP, visiting the providers/customers occurs during a tour. Because of ignoring tours in network design, location-allocation models in these systems will lead to higher costs [6]. Obviously, LRP is related to the classic location problem and VRP, both of which are expressed as special cases of LRP. If direct communication of all customers with a warehouse is necessary, then the LRP will be converted to a standard problem of location. On the other hand, if the position of warehouses is considered as fixed, LRP reduces to the problem of routing a simple vehicle.

In general, the LRP seeks to determine the location of warehouses and routes of vehicles for providing service to customers. Constraints include the capacity of warehouses, capacity of vehicles, route length, etc., to satisfy the demands of all customers and to minimize the total cost on the network, including routing costs, fixed costs of vehicles, fixed costs of facilities, and operational costs of facilities [7].

The first steps in LRP date back to the mid-1960s, when problems such as locating the shortest route arose. In the model presented by Hassanzadeh et al. [8], the costs of transportation were considered as part of the sum of the supply points from other points. The distances were weighted based on the volume of transportation between the points. Although the model has taken into account the transportation costs, it has not considered the constraints and complexities of the LRP model, which was far from a combined problem because, in the problem, the transportation cycles were ignored.

First real models of LRP were proposed in the late 1970s and early 1980s by several researchers like Jacobson and Madsen [9], [10], Or and Pierskalla [11], Nambiar et al. [12], Harrison [13], Laporte and Nobert [14], and Madsen [15].

During the recent years, researchers and logistics specialists have investigated more applied variants of LRP as well as more diverse methods of solving the problem. In general, progress on these articles coincided with the introduction of logistics and international trade growth, which were strongly dependent on the performance of the distribution system. Therefore, such problems are further applied and do not consider only academic aspects. A 2007 survey by Nagy and Salhi [16] highlighted that numerous researchers have explored the LRP by incorporating constraints related to warehouse or vehicle capacities. However, specific scenarios exist where vehicles or warehouses are assumed to have unlimited capacities. For instance, when the triangle inequality holds, vehicle routes must be reconstructed for each open warehouse. These variations typically frame the VRP as a node routing issue, where customers are treated as a subset of network nodes. Additionally, some studies, such as those by Ghiani and Laporte [17], Labadi [18], and Doulabi and Seifi [19], have examined the arc routing version of the problem.

Exact algorithms for the Capacitated Location-Routing Problem (CLRP) have emerged relatively recently. Belenguer et al. [20] proposed a branch-and-cut algorithm based on a linear binary model, enhanced by a new set of valid inequalities. Another approach, introduced by Akca et al. [21], utilizes the set partitioning framework and employs a column generation technique. Baldacci et al. [22] developed a method that decomposes LRP into multiple VRPs tied to a limited number of warehouses, leveraging dynamic programming and dual ascent strategies. Contardo et al. [23] later introduced a branch, cut, and price algorithm. These methods can efficiently solve problems with up to 50 customers and 5 to 10 warehouses, with the latter two also capable of handling larger instances involving 200 customers and 10 to 14 warehouses.

Glicksman and Penn [24] explored group LRPs, where warehouses and vehicles have no capacity constraints, and customers are grouped. Their approach selects a subset of open warehouses and tours covering all customer groups to minimize tour costs and fixed warehouse costs. In generalized VRP, each customer group has its own demand, and vehicles have limited capacities, with each group associated with a single warehouse. Harks et al. [25] analyzed worst-case performance ratios for various LRPs, such as uncapacitated depot LRP, cross-docking LRP, and prize-collecting LRP, by combining algorithms and reducing barriers across facilities. They achieved an approximation ratio of 4.38 for uncapacitated warehouse LRP. Carnes and Shmoys [26] studied K-LRP, where K uncapacitated warehouses with no fixed costs are considered in an open metric space. They developed a primal-dual scheme using Lagrangian relaxation to derive 2-approximation algorithms. Albareda-Sambola et al. [27] proposed a Tabu Search (TS) heuristic, starting with a solution derived from relaxed LP and rounding. The TS alternates between intensification at the routing level and diversification at the location level. Their results were compared using the AS benchmark set.

Derbel et al. [28] introduced an Iterative Local Search (ILS), which iteratively improves solutions through local search and perturbation operators. Their ILS reduced the average gap to lower bounds from 10.46% to 9.83% compared to Albareda-Sambola’s TS. In previous study [29], they extended ILS by embedding it into a GA, creating a memetic algorithm (GA-ILS) that slightly improved performance (9.81% vs. 9.83%).

Jarboui et al. [30] analyzed Variable Neighborhood Search (VNS) heuristics, encoding solutions as lists of customer routes per depot. Their best VNS combination marginally outperformed GA-ILS (9.70% vs. 9.81%) and was significantly faster. Prins et al. [31] developed the first matheuristic for CLRP, the Lagrangean Relaxation-Granular Tabu Search (LRGTS), separating location and routing decisions. This approach was later adapted by Ting and Chen [32] and Escobar et al. [33]. Duhamel et al. [34] replaced GRASP’s local search with Evolutionary Local Search (ELS), creating GRASP_ELS, which was further applied to CLRP and heterogeneous fixed fleet VRP by Duhamel et al. [35]. Yu et al. [36] implemented a simulated annealing heuristic for CLRP, encoding solutions as lists of warehouse sub-lists with fictitious zeros separating routes. Random moves included relocations, exchanges, and 2-opt moves. Jokar and Sahraeian [37] proposed another SA heuristic, starting with a greedy solution and improving it through local search. Jabal-Ameli et al. [38] introduced a Variable Neighborhood Descent (VND) algorithm for CLRP.

Ting and Chen [32] developed a Multiple Ant Colony Optimization (MACO) algorithm, decomposing CLRP into three decision levels: warehouse location, customer allocation, and capacitated VRP for each warehouse. Zarandi et al. [39] studied CLRP with fuzzy travel times and customer time windows, solving the model using Liu’s credibility theory and a simulated annealing algorithm. They tested the approach on a 20-customer, 5-warehouse example. Zarandi et al. [40] extended this work to include fuzzy customer demands and travel times, proposing a fuzzy chance-constrained programming model and an evolutionary method based on fuzzy c-means clustering and the sweep method. Their approach effectively solved problems with up to 100 nodes.

Mehrjerdi and Nadizadeh [41] addressed LRP with capacitated facilities and vehicles under uncertain demands modeled as triangular fuzzy variables. They used an evolutionary clustering algorithm and stochastic simulation to determine real customer demands, with improvements driven by an ant colony system. Their numerical experiments examined the risk of route failure and compared results with a commercial solver. Nadizadeh and Nasab [42] extended this work to include multiple vehicle routes. Golozari et al. [43] introduced an LRP variant with maximum travel time constraints, where customer demands, travel times, and service times were modeled as fuzzy variables. They used a fuzzy ranking function to linearize the model for small instances and a simulated annealing algorithm with a mutation operator for larger cases.

A traditional LRP typically assumes that customers only have delivery demands, focusing on how goods from open warehouses are distributed to customers using a fleet of vehicles. However, in real-world scenarios, customers often have both pickup and delivery demands, which may need to be fulfilled simultaneously. To address this, a specific variant of LRP has been introduced, known as the LRP with Simultaneous Pickup and Delivery (LRPSPD) [1]. This model has been applied in various reverse logistics contexts. Beverage Industry: LRPSPD is particularly relevant in the beverage sector, where companies are responsible not only for delivering products like beer and juice but also for collecting empty bottles for recycling and reuse [1]. Food Chain Stores: In the design of distribution systems for food chain stores, LRPSPD is utilized to manage the transportation of goods on pallets from distribution centers to stores, as well as the return of empty pallets to the distribution centers [1]. LRPSPD combines two core challenges: facility location decisions and vehicle routing with simultaneous pickup and delivery. As discussed in the work of Ismail Karaoglan et al., the flow-based model has proven effective in solving various LRPs. Consequently, this article adopts the flow-based approach to examine LRPSPD. In the flow-based model, additional variables are defined on the arcs of the graph, whereas the node-based model defines these variables based on the nodes of the graph.

We start with basic definitions of some well-known fuzzy numbers [8].

Definition 1: If a triangular fuzzy number ($LR$) with $\tilde{A}$ and membership function is $\mu_{\tilde{A}}$(x), the following equation holds:

\[\mu_{\tilde{\mathrm{A}}}(\mathrm{x})= \begin{cases}L\left[\frac{m-x}{a}\right] & x \leq m \\ R\left[\frac{x-m}{b}\right] & x \geq m\end{cases}\]

where, in the fuzzy number of $\tilde{A}=(m, a, b)_{L R}$, a denotes the left-hand side width and b denotes the right width (Figure 2).

Definition 2: Membership function of the trapezoid fuzzy number of $\tilde{a}=\left(a_1, a_2, a_3, a_4\right)$ is as follows:

\[\mu_{\tilde{a}}(x)=\left\{\begin{array}{cc}0 & x \leq a_1 \\\frac{x-a_1}{a_2-a_1} & a_1 \leq x \leq a_2 \\1 & a_2 \leq x \leq a_3 \\\frac{a_4-x}{a_4-a_3}& a_3 \leq x \leq a_4 \\0 & a_4 \leq x .\end{array}\right.\]

As shown in Figure 3, a triangular fuzzy number is a special case of a trapezoidal fuzzy number, where, $a_2=a_3$.

Definition 3: $\alpha$-cut and strong $\alpha$-cut for the fuzzy set of $\tilde{A}$ are shown by $\tilde{\mathrm{A}}_\alpha$ and $\tilde{\mathrm{A}}_\alpha^{+}$, respectively, with the following definition for $\alpha \in[0,1]$.

\[\begin{aligned}\tilde{\mathrm{A}}_\alpha & =\left\{x \mid \mu_{\tilde{\mathrm{A}}}(x) \geq \alpha, x \in \mathrm{X}\right\} \\\tilde{\mathrm{A}}_\alpha^{+} & =\left\{x \mid \mu_{\tilde{\mathrm{A}}}(x)>\alpha, x \in \mathrm{X}\right\}\end{aligned}\]
where, X is a general set.

Calculating $\alpha$-cuts for fuzzy numbers: $\alpha$-cut for a fuzzy number $LR$. For a fuzzy number $LR$ with input functions of $L$ and $R$, $\alpha$-cut is defined as follows (Figure 4).

\[\begin{aligned}& \alpha=L\left[\frac{m-x}{a}\right] \rightarrow \quad \frac{m-x}{a}=L^{-1}(\alpha) \rightarrow \tilde{\mathrm{A}}_\alpha^L=x=m-a L^{-1}(\alpha), \\& \alpha=R\left[\frac{x-m}{b}\right] \rightarrow \quad \frac{x-m}{b}=R^{-1}(\alpha) \rightarrow \tilde{\mathrm{A}}_\alpha^R=x=m+b R^{-1}(\alpha),\end{aligned}\]

$\alpha$-cut to a trapezoidal fuzzy number:

Let $\tilde{a}=\left(a_1, a_2, a_3, a_4\right)$ is a trapezoidal fuzzy number, $\alpha$-cut for $\tilde{a}$ is defined as follows:

\[\begin{aligned}& \alpha=\frac{x-a_1}{a_2-a_1} \Rightarrow \tilde{a}_\alpha^L=x=\left(a_2-a_1\right) \alpha+a_1, \\& \alpha=\frac{a_4-x}{a_4-a_3} \Rightarrow \tilde{a}_\alpha^R=x=a_4-\left(a_4-a_3\right) \alpha,\end{aligned}\]

Fuzzy numbers sum:

Let $\tilde{a}=\left(a_1, a_2, a_3, a_4\right)$ and $\tilde{b}=\left(b_1, b_2, b_3, b_4\right)$ are trapezoidal fuzzy numbers, then the following equation holds:

\[\tilde{a}+\tilde{b}=\left(a_1+b_1, a_2+b_2, a_3+b_3, a_4+b_4\right)\]

Comparison of fuzzy numbers:

The distance between two fuzzy numbers is shown by $D_{p, q}$, where index $p$ represents the analytical properties of the number (1$<$$p$$<$$\infty$), and $q$ is mental weight characteristic assigned to covered endpoints. If there is no reason for distinction between the sides of the fuzzy number, then $D_{p, 0.5}$ is recommended. Letting $q$ close to 1 means the right-hand side of fuzzy numbers is more favorable. Since importance of the endpoints covered by fuzzy numbers is assumed to be identical, the $q$ value is considered to be 0.5 (0$<$$q$$<$1).

\[D_{p, q}(\tilde{a}, \tilde{b})=\left[(1-q) \sum_{i=1}^n\left|a_{\alpha_i}^{-}-b_{\alpha_i}^{-}\right|^p+q \sum_{i=1}^n\left|a_{\alpha_i}^{+}-b_{\alpha_i}^{+}\right|^p\right]^{\frac{1}{p}} \]

If for two fuzzy numbers $\tilde{a}$ and $\tilde{b}$, there are n distinct numbers between zero and one for their $\alpha$-cuts, then the distance between the two triangular fuzzy numbers can be obtained from the following equation:

\[D_{2, \frac{1}{2}}(\widetilde{a}, \widetilde{b})=\sqrt{\frac{1}{2} \sum_{i=1}^n\left(a_{\alpha_i}^{-}-b_{\alpha_i}^{-}\right)^2+\frac{1}{2} \sum_{i=1}^n\left(a_{\alpha_i}^{+}-b_{\alpha_i}^{+}\right)^2} \]

To compare two fuzzy numbers $\tilde{a}$ and $\tilde{b}$ using $\alpha$-cuts, when they have positive values, we compare them with the fuzzy number of $MV=(0,0, \ldots, 0)$. In fact, we use $D(\tilde{a}, 0)$ and $D(0, \tilde{b})$ for comparison of the two numbers so that the number with a shorter distance with $MV$ is considered the smaller number.

Let $G=(N, A)$ be a fully directed network and $N=No{\cup}Nc$ be a set of nodes. No and Nc represent nodes associating to the potential warehouses and customers, respectively, and $A=\{(i, j): i, j \in N\}$ represents the set of arcs. Each arc $(i, j) \in N$ has a non-negative cost (distance) of $c_{i j}$, which holds in the triangle inequality, i.e. $c_{i j}+c_{j k} \leq c_{i k}$. Each potential warehouse $\mathrm{k} \in N o$ has the capacity of $CD_k$ and the fixed cost of construction is $FD_k$. An unlimited fleet of heterogeneous vehicles with a capacity of $CV_v$ and a fixed operating cost of $FV_v$ (including the cost of obtaining vehicles used on the road), is available to service customers. It should be noted that the assumption of unlimited dimensions of a fleet of vehicles with known capacities is a strategic point of view. Each customer $i \in Nc$ has a receipt demand $d_i$ and a pickup demand $p_i$, where $0 \leq d_i$ and $p_i \leq C V$. We also incorporated the following parameters in assumptions of the problem:

$l_v$: cost of the vehicle type of $v$ per each unit of distance

V: set of vehicles

In fact, the problem aims at determining the location of warehouses, allocating customers to opened warehouses, determining the route of related vehicles and allocating ideal vehicles to the relevant routes, and minimizing the total cost under the following constraints:

• Each vehicle uses at most one path.

• Each customer uses directly serves by a vehicle.

• All customers to a rout should only use a vehicle.

• Each path begins from a warehouse and ends to the same warehouse.

• The number of customers allocated to each open warehouse must not exceed the capacity of the warehouse so that for each open warehouse only one route is being used.

• Total load of vehicle at all points of the path does not exceed the capacity of the vehicle.

• Total sent demand and total delivered demand to customers allocated to each of the warehouses do not exceed the capacity of the warehouse.

Variables:

$x_{i j}=$ If a vehicle travel passes the direct route between $i$ and $j$, its value is 1, and otherwise is $0(\forall i, j \in N)$.

$y_k=$ If the warehouse $k$ is used, then its value is 1, and otherwise is $0 (\forall k \in No)$.

$z_{i k}=$ If customer $i$ is allocated to warehouse $k$, then its value is 1, and otherwise is $0 (\forall i \in Nc, \forall k \in No)$.

$f_{i v}=$ If the vehicle type $v$ is allocated to node $i$, then its value is 1, and otherwise is $0 (\forall i \in N, \forall v \in V)$.

Additional variables:

$D_{i j}$: The demand that must be delivered to customers on the route from node $i$ of goods that must be transported to the $\operatorname{arc}$ ($i, j$) (if a vehicle directly goes from node in to node $j, \forall i, i \in N$, then its value is 1 and otherwise is 0).

$P_{i j}:$ The sent demands of customers that are collected in route until the node ${i_w}$, which also incorporates the sent demand for the node $i_w$ in the arc ($i,j$) (if a vehicle directly travels in node i to node $j, \forall i, j \in N$, then its value is 1 and otherwise is 0).

Mathematical model:

Mathematical model of proposed problem is given below.

Model description:

Objective function (1) minimizes the total system costs, including shipping costs, fixed costs of warehouse construction, and using vehicles.

Constraints 2 and 3: These constraints are known as degree constraints. Constraint 2 ensures that each of the customers is met exactly once. Constraint 3 ensures that the numbers of input and output arcs to each node are equal. Constraint 4: This constraint shows that each customer is assigned only to one warehouse. Constraints 5, 6 and 7: These constraints prevent the creation of illegal routes, i.e., the routes that start from a warehouse and end at the same warehouse. Constraints 8 and 9: These constraints guarantee that total reception and delivery (separately) in each warehouse do not exceed the warehouse capacity. Constraint 10: This constraint shows that each customer is assigned only to 1 vehicle. Constraint 11: This constraint shows that nodes connected together take service only from 1 vehicle so that if node i is associated with the vehicle type v and i is connected to j directly, then j cannot belong to any other type of vehicle. Constraints 12 and 13: These constraints are flow inequalities for pickup and receive demand, and both of them prevent the creation of sub-tours and ensure that pickup and receive requests of each of the customers are met. Constraint 14: This constraint emphasizes the concept that the total load in each arc should not exceed the capacity of the vehicle. Constraint 15: This constraint ensures that the entire shipped load from each open warehouse is exactly equal to the delivery demand of all customers assigned to the warehouse. Constraint 16: This constraint ensures that the entire load shipped back to the open warehouses must be equal to zero. Constraint 17: This constraint ensures that the total load delivered to each open warehouse is exactly equal to the sent demand of customers assigned to the relevant warehouse. Constraint 18: This constraint ensures that the total sent load in the vehicle when every dispatch from the open warehouse should be zero. Constraints 19, 20, 21 and 22: These are bounded constraints for decision variables. Relations (23) to (26) are the sign and type of variables.

As we mentioned earlier, parameters relating to the customer's request, i.e., di and pi, are considered to be fuzzy. In addition, because of equality constraints, i.e., constraints 12, 13, 15 and 17, inevitably to solve the problem, decision variables $D_{ij}$ and $P_{ij}$ should be transformed into fuzzy. The fuzzy mathematical counterpart model of the proposed problem is given below.

GA can be considered as a directional stochastic optimization method that gradually moves toward an optimal point. Before a GA can be implemented, the first chromosome structure or an appropriate presentation should be found for a related problem. Also, a fitness function should be devised to assign a value to each encoded solution. During the implementation process, the parents will be selected for reproduction, and through crossover and mutation operators, new offspring will be produced. This process will be repeated several times to produce the next generation of the population. Then the population will be investigated; if convergence criteria are met, the process will be terminated.

Chromosome structure:

$A((Nc+No)*(No+No))$ square matrix will be created from the path between nodes ($xij$) as follows: The rectangular sub-matrix $(Nc*(Nc+No)$ will be created, so that each row of it has only one $cij$ equal to 1, and columns do not have more than one xij equal to 1. Finally, for the remaining rows related to distribution centers ($No$), the cells of the matrix ($xij$) will be filled such that correct paths are created. Therefore, the chromosomes of our GA are a rectangular matrix $(Nc*(Nc+No))$, where each row of the matrix (from 1 to $Nc$) is considered a gene algorithm.

Example 1:

3 customers and 2 distribution locations are given. Numbers 1 through 3 are dedicated to the customer, and numbers 4 and 5 are dedicated to the location of the distribution centers, and a square matrix will be created as follows. Since each cell of the matrix represents one edge of the path, then some of the table cells cannot be equal to 1; these cells are indicated with a black color. The cells located on the main diameter cannot be equal to 1 because no travel can happen from one place to the same place. Also, the square cells located at the end of the matrix cannot be equal to 1 because there is no path between distribution centers.

As can be seen in Figure 5, first $(Nc*(Nc+ No))$ submatrix values and then the values related to the distribution centers will be allocated. Therefore, in this question, our chromosome algorithm will be $(Nc*(Nc+ No))$ matrix that is illustrated in Figure 6.

Initial population: In order to make the initial population, one cell will be randomly selected for each row, and a value equal to 1 will be assigned to it. In a row, which one of the cells “is qualified to receive a value equal to 1”? The cells are also qualified to receive a value equal to 1, so that the total of the values related to that column becomes zero.

Example 2:

In Figure 7, all cells qualified to receive a value equal to 1 for row 3 are shown below.

Selection operator: “Roulette wheel" Operator was used in order to select the parents for the next generation. In this cycle, more qualified chromosomes had more chance to be selected for reproduction.

Crossover operator: In a cross, depending on the kind of crossover operator, two or more points located between rows of the matrix will be selected, and then the rows of the two matrixes will be replaced. In each column, more than one cell cannot be equal to 1. Therefore, if after crossover, summing each column became more than 1, the rows that have caused this violation should be regarded equal to zero, and in each one of these rows, one of the cells that is qualified to be equal to 1, will receive a 1 value.

As it can be seen in Figure 8, for the first child, in the second row, more than one cell is equal to 1. Therefore, in the crossed row, the qualified cell is to be 1 and is marked with a red circle, will receive a 1 value, and the matrix will be modified.

Mutation: For mutation, one or more rows of the $(Nc*(Nc+ No))$ matrix will be selected. Then a value of zero will be assigned to all cells of these selected rows. Then, a qualified cell in each selected row is randomly chosen, and its value is set to one.

Note: The first line to which the mutation operation is applied will be selected randomly. For Example 1, if rows 1 and 3 of a matrix are going to suffer a mutation, the first row will be selected randomly, and the mutation will be applied to it. In the mutation process, for the first row to which the mutation is applied, the cell that previously had a value of 1 will not be considered qualified to become 1 again. However, this rule does not apply to other rows. Subsequently, the qualified cells are illustrated in Figure 9, for example.

The fitness function is similar to the objective function, with the addition of M times the number of unmet constraints, where M is a very large number. The goal is to minimize this value.

The GA used to solve the non-fuzzy model is also applied to the fuzzy model. However, for the calculation of the GA's fitness function, the operators described in the fuzzy section are applied, including the summation of fuzzy numbers and the comparison of two fuzzy numbers. The crucial aspect of fuzzy computing in the GA lies in determining whether the constraints are satisfied. By using the comparison of fuzzy number's formula derived from the paper of Hassanzadeh et al. [8], the problem constraints are investigated. Constraints that are satisfied receive a zero penalty in the fitness function, while unsatisfied constraints impose a penalty to guide the algorithm toward a more feasible solution space. Additionally, to calculate the distance in urban areas, the norm-1 between two points is applied, as used in this model as well.

Based on the values of $x_{i j}$, the route variables obtained indicate which customers are assigned to which warehouses (denoted as $z_{i k}$). Consequently, it can be concluded that only the warehouses assigned to at least one customer are established, while those with no assigned customers cannot be established (denoted as $y_k$).

Method of Obtaining the $f_{i \mathrm{k}}$ Variable: This variable indicates which customer is assigned to which vehicle, given that only one vehicle can be assigned to each route. Initially, all customers on the same route-defined as those traveling from the moment a certain warehouse is departed until the moment they return to that warehouse-are assigned to the same vehicle. If we arrange all types of the vehicles from small to large, the first type of vehicle in which the capacity is greater than the sum of $\widetilde{p}_l \mathrm{s}$ and greater than the sum of $\widetilde{d}_l \mathrm{s}$ will be assigned to that route.

Then, based on the values obtained for the variables $\widetilde{D_{l j}}$ and $\widetilde{P_{l j}}$ from constraints $11,12,14,15,16$, and 17, the vehicle for the route is selected to satisfy constraint 13. If no vehicle meets this constraint, a vehicle with a capacity one unit greater than the previously selected one is chosen. This process continues until a vehicle satisfies all the constraints. However, if no vehicle meets the constraints, even after considering the largest available vehicle, the route is deemed invalid, and the route matrix must be reviewed to re-select a new route.

The following terms are used in the subsequent discussion:

$x_c$: horizontal coordinate (abscissa) associated to each customer,

$y_c$: vertical coordinate associated to each customer,

$x_o$: horizontal coordinate (abscissa) associated to each customer warehouse,

$y_o$: vertical coordinate associated to each warehouse,

$C D$: the capacity of each warehouse,

$FD$: fixed cost of constructing every warehouse,

$d l$: lower bound of the fuzzy number of the times that each customer has received,

$d$: the fuzzy number of the times that each customer has received,

$d u$: the fuzzy upper limit number of the times that each customer has received,

$p l$: the fuzzy lower bound of the number deliveries for each customer,

$p$: the fuzzy number of the number deliveries for each customer,

$pu$: the fuzzy upper bound of the number deliveries for each customer,

$v$: index of relevant vehicle,

$C V$: the capacity of each vehicle,

$F V$: the fixed cost relating to each vehicle,

$l$: the cost of each vehicle per unit of distance.

Note: The indices whose values exceed the number of customers correspond to warehouses. For example, in Example 3, $x_{19}=1$ means that the route from customer 1 to warehouse 3, i.e. $3=9-6$, is passes by one vehicle.

Example 3:

6 customers and 3 potential warehouses.

The problem's parameters and initial data are presented in Table 1, while the problem's diagram is shown in Figure 10. In the diagram, potential warehouse locations are represented by squares, and customer locations are depicted by circles.

Answer of Example 3 on non-fuzzy status:

\[\begin{aligned} & \mathrm{W}^*=832 ; x_{94}=x_{41}=x_{19}=1, \quad f_{91}=f_{41}=f_{11}=1; \quad x_{83}=x_{38}=1, \quad f_{81}=f_{31}=1;\quad x_{76}= \\ & x_{65}=x_{52}=x_{27}=1, \quad f_{72}=f_{62}=f_{52}=f_{22}=1\end{aligned}\]

Answer of Example 3 on fuzzy status:

\[\begin{aligned} & \mathrm{W}^*=1721; \quad x_{95}=x_{52}=x_{29}=1, \quad f_{94}=f_{54}=f_{24}=1 ; x_{84}=x_{43}=x_{38}=1, \quad f_{83}=f_{43}=f_{33}= \\ & 1 ; \quad x_{71}=x_{16}=x_{67}=1, \quad f_{73}=f_{13}=f_{63}=1\end{aligned}\]

As shown in Figure 11 and Figure 12, the networks in the fuzzy and non-fuzzy states are completely different.

In this study, the LRPSPD model introduced by Karaoglan et al. [1] was investigated. The primary objective was to improve the model's applicability by relaxing the assumption of a homogeneous vehicle fleet. Real-world applications typically involve a heterogeneous fleet, thus this variation was incorporated into the model. This modification resulted in the introduction of different vehicle capacities and associated costs, which necessitated the inclusion of new constraints and slight adjustments to the objective function. A further critical aspect addressed in the study was the uncertainty in customer demands. In practical scenarios, detailed information regarding customer delivery and pickup demands is frequently unavailable. However, expert estimates, based on prior experience, are often relied upon. To handle this uncertainty, customer delivery demands ($d_i$) and pickup demands ($p_i$), such as the return of goods to distribution centers, were treated as fuzzy variables. This approach significantly enhanced the model’s realism, making it more representative of real-world conditions. The model was solved using GA, and the results indicated that the total cost in the fuzzy model was higher than in the non-fuzzy model. Despite this, the fuzzy model ensured that all constraints were satisfied, demonstrating its practical advantage. The results also highlighted that ideal planning outcomes, based on certainty, are often unachievable due to the inherent uncertainties in real-world data. By incorporating fuzzy logic, the uncertainty in the planning process was effectively addressed, which resulted in a more robust model aligned with real-world administrative practices. The findings underscore the importance of accounting for uncertainty in LRP. The inclusion of fuzzy logic not only improved the model’s practical relevance but also provided a more accurate representation of complex, uncertain environments. This work contributes to the field by extending traditional deterministic models, offering a more flexible and adaptable approach for solving LRPs in real-world logistics networks.