1. INTRODUCTION
Transportation is the pillar of a communication bridge passed by various community sections to move toward sustainable development. In today’s world, transportation has become a fundamental part of the economy and has affected the process of economic development. In addition, it is the basis for trade exchanges and the key to economic and social development. In fact, transportation acts as a vein that connects another service, commercial, industrial and agricultural components to each other at the national and international level, exerting deep impacts on the growth and development of these sectors, the performance of each of which plays a direct role in the progress and regression of a country's economy. Some of the most important causes of increased transportation costs include crowded roads (Sankaran and Wood, 2007), car accidents, travel delays (Steimetz, 2008) and energy consumption. Given the effect of distance on all of the mentioned areas, it could be expressed that even the slightest improvement in the path can lead to a considerable decrease in costs (Santos et al., 2010). The importance of this issue has drawn the attention of many researchers in vehicle route design. Meanwhile, the tracking issue has been a consi derable part of the process. In addition to determining the route of the service vehicles, deciding about the location of depots (the start and ends of tours) is significantly important in designing distribution systems. Over the past few years, making effective, reliable, and flexible decisions about these two topics has become a major concern to managers (Nadizadeh et al., 2011).
Different researchers have often criticized the separate assessment of these two key areas. For instance, Salhi and Rand (1989) showed that solving the locationrouting problem (LRP) while overlooking the routes could lead to a suboptimal solution. To solve this issue, some examples were proposed, where decisions related to locationrouting problems were made simultaneously (Derbel et al., 2012;Karaoglan and Altiparmak, 2015). Through simultaneous assessment of the facility location problem (FLP) and vehicle routing problem (VRP) under limitations such as the capacity limit of vehicles and facilities, a solution is proposed by which the demands of customers are responded to with the minimum overall cost (including the transportation, depo construction, and vehicle fixed cost). The LRP is used in various areas, including newspaper delivery (Jacobsen and Madsen, 1980;Madsen, 1983), and the distribution systems of beverage products (WatsonGandy and Dohrn, 1973). In an LRP, customers only demand delivery, and the goal is to determine how to distribute the items using vehicles based on active depots. One of the hypotheses that can make the problem more applicable to some realworld problems is that pickup be demanded by customers in addition to delivery, and these two types of demands be responded to simultaneously and by one vehicle. To solve this issue, Karaoglan et al. (2011) defined problems known as location routing problem with simultaneous pickup and delivery (LRPSPD) as a generalized form of LRP and VRPSPD (vehicle routing problem with simultaneous pickup and delivery) (Karaoglan et al., 2011). A brief explanation of the problem is simultaneous determining the location of active depots (with predefined capacity and vehicle fixed cost) and the route of service vehicles (with predefined capacity and vehicle fixed cost) to simultaneously satisfy the pickup and the delivery demands of each customer at the lowest overall cost. This method is used in distribution systems of beverage products, where vehicles should pick up empty bottles of customers in addition to the delivery of products to these individuals.
Since the LRP and LRPSPD share the same feature of being an NPhard problem (Vincent and Lin, 2014), precise optimization methods lack the efficiency to solve medium and largescale instances. Frequently showing their ability to achieve the optimal (nearoptimal) solution of NPhard problems at a faster pace, metaheuristic algorithms have been recognized as one of the best options available to solve such problems. The current study aimed to present a new approach to solving LRPSPD using the wellknown particle swarm algorithm (PSO). The remainder of this paper is structured as follows: section 2 reviews the literature and section 3 defines the problem and presents a mixedinteger programming (MIP) model. In addition, section 4 describes the approach proposed to solve the problem and section 6 evaluates the performance of the proposed algorithm using a numerical example. Finally, section 6 concludes.
The LRPSPD was first introduced by Karaoglan et al. (2011) who proposed a problemsolving method after presenting a mathematical model, which encompasses five different decision variables, using branch and cut and simulated annealing (SA) optimization techniques, as well as several valid inequalities. Using the samples extracted from the literature to assess the performance of the new algorithm, they concluded that some samples consisting of a maximum of 88 customers and 8 potential depots could be solved in a reasonable time. Later on, Karaoglan et al. (2012) presented two polynomialsize mixedinteger linear programming formulations (one is nodebased and the other one is flowbased) and a family of valid inequalities to strengthen the formulations. They also proposed a twophase heuristic approach based on SA, tp_SA to solve the largesize LRPSPD, and two initialization heuristics to generate an initial solution for the tp_SA. According to the results, the flowbased formulation performed better than the nodebased formulation in terms of the solution quality and the computation time on smallsize problems (1030 customers). However, the nodebased formulation yielded competitive lower bounds in a reasonable amount of time on mediumsize problems (50100 customers). Moreover, the SA algorithm reached a mean distance of 1% and a maximum distance of 6% with the best lower limit in mediumsized samples.
Vincent and Lin (2014) proposed a multistart simulated annealing (MSA) algorithm for solving LRPSPD, which incorporated a multistart hillclimbing strategy into a simulated annealing framework. They tested the MSA algorithm on 360 benchmark instances to verify its performance, and the results were indicative of a significant enhancement in the performance of the traditional singlestart SA algorithm by using the multistart strategy. The aforementioned scholars carried out another research later on, where they proposed an SA heuristic for the LRPSPD by using the power of the SA algorithm (Vincent and Lin, 2016). The algorithm’s performance was assessed using eight sets of samples. According to the results, the proposed SA effectively solved LRPSPD. GhatrehSamani and HosseiniMotlagh presented a mixedinteger linear programming model for a twoechelon LRPSPD (Ghatreh Samani and HosseiniMotlagh, 2017). In the mentioned research, uncertain demands were taken into account and a fuzzy programming approach was exploited to deal with fuzzy parameters. The researchers devised a combined heuristic method based on simulated annealing (SA) algorithm and genetic algorithm (GA) for solving the presented model, the performance of which was assessed using different numerical examples. Wang and Li (2017) studied the decrease of carbon emissions for LRP with a heterogeneous fleet, simultaneous pickupdelivery, and time windows. In addition to presenting a MIP model for the problem, they designed a twophased hybrid heuristic algorithm to generate an initial solution for variable neighborhood search (VNS). These scholars also exploited the genetic algorithm and improved the optimization power of the algorithm by incorporating the idea of the SA algorithm into the VNS framework.
Nadizadeh and Kafash (2017) modeled the fuzzy capacitated locationrouting problem with simultaneous pickup and delivery demands (FCLRPSPD) based on the fuzzy credibility theory. In addition, they developed a greedy clustering method (GCM) comprising of four iterative phases to solve the FCLRPSPD. After clustering the customers, the center points of the clusters for selecting the suitable depot(s) and allocating clusters to the selected depot(s) were determined in the final phase by the ant colony algorithm to create a path between depot(s) and the allotted clusters. In another research, Nadizadeh and Hosseini Nasab (2019) used the mentioned greedy clustering method (GCM) to solve the problem after presenting a nodebased MIP formulation of the CLRPSPD (capacitated locationrouting problem with simultaneous pickup and delivery). In the latest research, Martins et al. (2002) proposed a mathematical and a heuristic for solving LRPSPD. The mathematical model was employed to reduce the last delivery time. In addition, fleet capacity constraints were considered, and a savingsbased approach was exploited to solve the problem. Results showed that the proposed heuristic could detect feasible and competitive solutions in a very short computational time. In a research, Nadizadeh and Kafash (2017) proposed a waiting strategy for the VRPSPD, for which they first presented a mathematical model where a great variety of products was delivered to customers by a distribution system, and their recycled materials were received. The researchers used a genetic algorithm to solve this mathematical model. The efficiency of the algorithm was shown by solving various numerical examples.
2. METHODOLOGY
The LRPSPD was defined on a complete directed graph (G) encompassing a set of nodes (N) (e.g., depot establishment nodes [N_{0}] and customer nodes [N_{c}]) and a series of arcs (A). A nonnegative number (C_{ij}) was attributed to each of the arcs (e.g., the I,j arc), which was the cost of travel from the ith node to the jth node (or the distance between the ith node and the jth node). The numbers should be chosen in a way that they are applied to a triangle inequality (i.e, (C_{ij} + C_{jk} ≥ C_{ik}) ). In addition, a capacity (CD_{k}) and fixed cost (FD_{k}) were considered for each depot establishment node. A fleet comprising of an unlimited number of identical vehicles with similar capacities (CV) and operation fixed costs (FV) was based in depots to serve the customers (each with a delivery [D_{i}] and pickup [P_{i}] demand). The purpose of locating depots is to allocate customers to the built depots and to determine the route of service vehicles with the lowest total cost considering the following hypotheses:

•Each device must be used in a maximum of one direction.

•Each customer must be served by only one vehicle.

•Each path should start and end with a single depot.

•A total load of a device should not exceed the capacity of that device at any point on the route.

•The total delivery demands of customers allocated to a depot should not exceed its capacity.

•The total pickup demands of customers allocated to a depot should not exceed its capacity.
Afterwards, we performed the flowbased MIP model (13), which includes five types of decision variables, as presented below:

Y_{k} is one if the depot is established at the kth location; otherwise, it is zero.

X_{ij} is one if a vehicle directly travels from the ith node to the jth node; otherwise, it is zero.

Z_{ik} is one if the ith customer is allocated to the kth depot; otherwise, it is zero.

U_{ij} equates the number of remaining delivery demands after leaving the ith node if the vehicle directly travels from the ith node to the jth node; otherwise, it is zero.

V_{ij}, equates the number of collected pickup demands up to the ith node if the vehicle directly travels from the ith node to the jth node; otherwise, it is zero.
In the mentioned model, formulation (1) minimized the total costs of transportation, depots and vehicles. Constraints (2) ensured that each customer was met exactly once and the set of constraints (3) ensured that the number of arcs entering each node was equal to the number of arcs leaving that node. It is notable that each customer had to be allocated to just one depot due to constraints (4). In addition, conditions 57 prevented the formation of unauthorized routes that did not start and end in the same depot. Moreover, constraints (8) and (9) ensured that the total pickup and delivery demand of customers allocated to a depot would not exceed the depot’s capacity. On the other hand, constraints (10) and (11) maintained the flow of pickup and delivery demands in addition to removing subtours. Based on condition (12), the total load passing through an arc must not exceed the vehicle capacity. Constraints (13) ensured that the total delivery load taken out of an active depot was equal to the sum of customer delivery requests allocated to that depot. Constraints (4) claimed that the total load returned to active depots must be equal to zero. In addition, constraints (15) were added to the model to ensure that the total input pickup load of an active depot was equal to the total pickup demands of customers allocated to that depot. Furthermore, constraints (16) ensured that the total pickup load removed from the active depots is zero. In addition, conditions 17 20 determined the lower and upper bounds of additional variables (U & V). Due to constraints 2123, the main variables (Y, X, and Z) can only take the values 0 and 1, and constraints (24) limited the additional variables to nonnegative values.
3. RESULTS AND DISCUSSION
The PSO is a global optimization method used to deal with problems whose answer is a point or surface in an ndimensional space. Some hypotheses are made in such a space and an initial velocity is assigned to the particles. In addition, communication channels are considered between particles, and the particles move in the response space. The results obtained are calculated in each timespan based on “competency criteria”. Over time, the particles accelerate toward those particles that have higher competency criteria and are in the same communication group (Sangaiah et al., 2020;Piotrowski et al., 2020;Michalewicz and Fogel, 2013). The main advantage of this method over other optimization strategies is that a large number of swarm particles makes the method flexible in the face of a local optimal response problem. Figure 1 depicts some examples of the particles’ movement pattern in the search space. In the picture, the primary location of the particles is the upper left corner of the image, which is twodimensional in the search space. The swarm will eventually converge with the repetition of the algorithm, which is observed in the lower corner on the right side of the image.
Each particle has a location that determines its coordinates in a multidimensional search space. Over time, the position of the particle changes with its movement. x_{i}(t) determines the position of the ith position in the tth time. Each particle needs a velocity to move in space. v_{i}(t) determines the velocity of the ith particle in the tth time. A new position can be considered for the particle by adding velocity to its location. The equation used to update the particle’s position is shown in (25):
The suitability of a particle’s position in the search space is assessed by a fitness function. Particles can remember the best situation they have ever been. The best individual performance of a particle or the best position met by the particle is called y_{i} (in some algorithms, y_{i} is also called pbest) and particles can know about the best position met by the entire group, which is known as ${\widehat{y}}_{i}$ (in some algorithms, ${\widehat{y}}_{i}$ is also called gbest). The particle velocity vector in the optimization process reflects the empirical knowledge of the particle and the information of the particle community.
Cognitive component: is the best solution obtained by a particle independently.
Social component: is the best solution detected by the entire group.
Piotrowski et al. (2020) showed a sample of flowchart of the PSO algorithm:
There are two main models for the standard PSO algorithm, the velocity vector of which is estimated based on the two cognitive and social components. The two models are known as the lbest PSO and gbest PSO, which are different in terms of neighborhood size considered for each particle. This is further explained below:
The gbest PSO model: in this algorithm, the neighbor is a particle of the whole group of particles, and the topology required to connect the particles to each other is similar to the star topology. In this algorithm, ${\widehat{y}}_{i}$ is the best position met by the entire group (88, 92, 93, 9599). The position of each particle is evaluated by the particle suitability criterion function. If the new position of the particle leads to a more optimal value, the position of the particle will change; otherwise, the position of the particle is constant. The value of ${\widehat{y}}_{i}$ is also equal to the most optimal value of the entire particles in each repetition. The particles are first scattered in the search space and then each particle converges to a point at a certain velocity. The velocity of a vector is Nddimensional, where v_{ij} exhibits the velocity of the jth element of the velocity vector of the ith particle. Therefore, the velocity of the ith particle is updated based on the formulation (26) (22).
Where W determines the inertia weight, C1, and C2 are constants. The term “inertia weight” was first introduced by Shi and Eberhart in 1998. This weight controls the impact of the previous velocity in calculating the new velocity. The higher the value, the higher the public search, and the lower the weight, the higher the local search.
Figure (2) shows how to update the particle velocity.
There should be a connection between the fixed values and the existing weight in order to ensure a solution. Otherwise, it causes particles to diverge due to a series of behaviors.
The position of the ith particle is updated by the formulation (28):
Particle position updating leads to the updating of the velocity and position of each particle by the mentioned formulations, and this process is repeated until reaching the maximum number of repetitions or having a nearzero updated velocity and estimating the performance of each particle by the competency criterion.
The lbest PSO model: in this model, a particle has only the ability to communicate with a number of particles in its neighborhood to update its velocity. In addition, N_{i} is the total number of neighboring particles of one particle? The ring topology is used to communicate between particles, and the best position met by neighbors of the particle is called ${\widehat{y}}_{i}$. Other neighbors are shown by l.
The formulation for updating the position of particles is similar to the previous method and does not change. The neighbors of a particle practically determine its social behaviors. Therefore, neighborhood topology is an important and useful topic.
In this research, the gbest PSO model is used to optimize the problem. In addition, Figure 4 shows the standard pseudocode of this algorithm:
Consider a problem with five depot nodes (nodes 15) and 20 customer nodes (nodes 625). In the present approach, a sample solution is exhibited based on Figure 1.
In Figure 1, the lines indicate the direction of movement of vehicles. The first number of each line (from the left side) shows the number of the depot node of the start and end of the travel, and the rest of the numbers determine the order of customer nodes traveled by a vehicle. According to the figure:

* The depot nodes 2 and 3 are the only active ones.

* Four vehicles are considered for serving the customers.
The trajectory of the first and second vehicles is 2 11 15 14 6 2 and 3 7 19 20 3, respectively.
This section analyzes the numerical results, and the main goal is to assess the GAMS’s ability to solve the mathematical model. Another objective is to determine the power and quality of the PSO algorithm in an approximate solution of the problem ahead. Initially, 10 samples are randomly generated in different dimensions. The dimensions of each problem are defined according to the limits of its indices. Therefore, the dimensions of each of the 10 generated examples are presented in Table 1, where l introduces customers and k shows the number of depots.
Information about each of the examples A1A10 is inserted in the GAMS software and in the MATLAB software as the input of the PSO algorithm and their output is recorded. It should be noted that each software has a time limit of 3600 seconds to solve each problem. Table 2 summarizes the results of solving the sample examples with the help of the mentioned software.
In Table 2, Z is the value of the objective function obtained from the desired solution method. T represents the software solution time and GAP is the percentage difference between the optimal value and the best value found by the PSO algorithm. Moreover, GAP is calculated as Equation 30.
Table 2 presents the results of the exact solution of the mathematical model and the solution using the PSO algorithm. Due to the time limit of 3600 seconds, GAMS software could not solve examples A8A10. On the other hand, regarding examples A6 and A5, the program is stopped and the best answer to the problem is provided after 3600 seconds, which is due to the high complexity of the mathematical model and problem, which inhibits its solving at larger scales with the help of the GAMS. Meanwhile, the average solution time in MATLAB software and PSO algorithm is less than 30 seconds and all examples are executed in less than one minute. In figures (35), the solution times of the two methods are compared in different examples.
As observed in Figure 6, the process of increasing the dimensions of the problem in GAMS software is exponential, which shows the Nphard nature of the problem under study. Now, careful examination of the solution time trend in the PSO algorithm shows that the slope of increasing the solution time is very low. Meanwhile, the average error of the PSO algorithm is 0.7%. In other words, the PSO algorithm imposes a 0.7% error resolution time on the decisionmaker with a very sharp reduction in time, which is very convenient compared to the reduction of the solution time. Figure (7) shows the error value of the PSO algorithm method on each example.
As shown in Figure 8, the maximum error generated by the PSO algorithm is 2.01%. On the other hand, the amount of trend error has not increased with increasing dimensions of the problem, which means that the PSO algorithm will provide the best possible answer with a reasonable error if there is a need for solving problems with much larger dimensions. As such, the PSO algorithm designed for this research can be introduced as a practical tool in solving LRPSPD.
4. CONCLUSION
The present study specifically focused on the optimization of location routing problem with simultaneous pickup and delivery, for which a new mathematical model was developed. In the model, a decision such as allocation of customers to depots, routing vehicles, and determining the volume of satisfied demands (of both pickup and delivery type) were optimized simultaneously. The PSO algorithm model was proposed to solve the model. In addition, several samples were implemented to assess the efficiency of the PSO algorithm, and its performance was compared in GAMS by accurate solving of the mathematical model. According to the results, the proposed algorithm could optimize the LRPSPD problem in a very short time. Meanwhile, the mean error of this algorithm was only 0.7%, which indicated the high efficiency of this solution method. In order to expand the research, it is recommended that uncertainty in important parameters of the mathematical model be considered, especially the level of received pickup and delivery demands. In addition, it is suggested that robust optimization methods and feasibility planning be used to deal with this uncertainty.