Reverse Logistics Location Based on Energy Consumption: Modeling and Multi-Objective Optimization Method

The low-carbon economy, as a major trend of global economic development, has been a widespread concern, which is a rare opportunity to realize the transformation of the economic way in China. The realization of a low-carbon economy requires improved resource utilization efficiency and reduced carbon emissions. The reasonable location of logistics nodes is of great significance in the optimization of a logistics network. This study formulates a double objective function optimization model of reverse logistics facility location considering the balance between the functional objectives of the carbon emissions and the benefits. A hybrid multi-objective optimization algorithm that combines a gravitation algorithm and a particle swarm optimization algorithm is proposed to solve this reverse logistics facility location model. The mobile phone recycling logistics network in Jilin Province is applied as the case study to verify the feasibility of the proposed reverse logistics facility location model and solution method. Analysis and discussion are conducted to monitor the robustness of the results. The results prove that this approach provides an effective tool to solve the multi-objective optimization problem of reverse logistics location.


Introduction
China is the world's most populous country, with a population of more than 1.41 billion as of 2021. At present, the academic circles have extensive discussions on the resource utilization of waste products, and the research on the reverse logistics of waste products is gradually being carried out with the continuous extension from both the academic world and industries in recent years [1,2]. According to the findings of Rogers and Tibben-Lembke, the total logistics cost amounted to USD 862 billion in 1997, and the total cost spent in reverse logistics is enormous and amounted to approximately USD 35 billion, which was around 4% of the total logistics cost in the same year [3]. The concerns about energy saving, green legislation and green manufacturing are increasing [4].
The concept of reverse logistics was put forward by Stock in the 1990s, which aims at high efficiency and environmental protection and which maximizes the application value of products by optimizing the operation structure [4]. With the ever-rising needs of reverse logistics, firms possessing optimal planning of return routes, inventory and warehouse layouts for returned products are more competitive than the others. Reverse logistics has changed the original resource flow mode, and the single logistics system has evolved into a closed-loop system because of the participation of reverse logistics. Reverse logistics, including product return, maintenance, remanufacturing, waste disposal and other processes, is the mirror reverse process of traditional logistics. The purpose of this method is to extract and utilize the value in the waste repeatedly and to improve the production chain [5]. Regarding distribution paths, the routing and scheduling of reverse logistics are more complex than those of forward logistics. The routing of reverse logistics starts from the designated regional distribution center to the centralized return center and then to manufacturers for remanufacturing or to manufacturers directly without passing through the centralized return center. Due to the uncertainty of return quantities, the physical flow channel is more complicated than that in forward logistics. Since different areas have various product return rates, the locations of collection points can significantly affect the efficiency of recycling [3].
A series of systemic theoretical explorations and literature reports for reverse logistics location have been investigated in recent years [5,6]. For example, Savaskan et al. established a decentralized decision system with three schemes and concluded that the agent nearest to the customer is the best one to undertake product recycling [7]. Park et al. analyzed the influences of different subjects on recycling price and profit and built a closed-loop supply chain model with a single recycling channel [8]. Liu et al. investigated the competition between formal recycling channels and the informal, concluding that high-quality recycled products were more attractive to the two recycling channels, while the non-normal recycling channels were willing to recycle at a higher price [9]. Tasbirul et al. estimated the economic value generated in the recycling logistics of electronic waste gas based on the sales inventory life model and provided suggestions for decision makers [10]. In addition, construction of a reverse logistics model is an important part of the research on reverse logistics location. Shih et al. studied the recycling model of waste household appliances with the aim of minimizing the recycling cost and analyzed the influence of recovery rate and storage state on the recycling network [11]. Lee et al. built a two-stage recycling model and carried out an entry study on the problem of uncertain recycling quantity [12]. To solve the multi-objective optimization problem of reverse logistics location, some evolutionary algorithms have been applied, e.g., the ant colony algorithm [13], simulated annealing algorithm [14], particle swarm optimization [15] and the genetic algorithm [16].
In summary, most of the existing research literature starts from cost recovery or benefit recovery and takes the maximum benefit or minimum cost as the objective function. However, economic benefit is only one of the benefits of reverse logistics, and how to use recycling behavior to produce considerable environmental benefits cannot be ignored. It is contrary to the original intention of the circular economy to take economic benefit as a single goal. Thus, it is necessary to take energy consumption or other indicators that can reflect the utilization of resources into consideration and to consider the balance between each indicator and the profit of each objective. In addition, how to efficiently achieve optimal design under multiple factors and objective functions has become a significant research topic and can be regarded as a typical multi-objective optimization problem. This paper formulates a double objective function optimization model of reverse logistics facility location considering the symmetry among the functional objectives of the carbon emissions and benefits. A hybrid multi-objective optimization algorithm that combines a gravitation algorithm and a particle swarm optimization algorithm is proposed to solve this reverse logistics facility location model. Subsequently, an empirical case of a mobile phone recycling logistics network in Jilin Province is applied to verify the feasibility of the proposed reverse logistics facility location model and the multi-objective optimization method. Analysis and discussion are conducted to monitor the robustness of the results.
The paper is organized as follows: Section 2 introduces a reverse logistics facility location model with the functional objectives of addressing carbon emissions and the resulting benefits. In Section 3, we develop a hybrid multi-objective optimization algorithm, combining the theories of the gravitation algorithm and the particle swarm optimization algorithm. An empirical case of the mobile phone recycling logistics network in Jilin Province is applied to verify the effectiveness of the model and solution methodology in Section 4. Analysis and discussion are conducted in Section 4. Finally, the conclusions of this paper are drawn in Section 5. To make the model easy to quantify and to make the optimization algorithm convenient to process, the following assumptions are made for the waste product recycling network: Assumption 1. Recyclers work with product sellers, so there is no need to build another collection point. Product sales points of each city are recycling points. This model is a simulation of waste product recycling in a province, so the recycling logistics in the city are not considered. That is, the logistics costs from the consumer to recycling points are ignored.

Assumption 2.
It is possible for a city to build both storage points and specialized treatment plants, and the recycling plant has the capacity to dispose of dismantled waste, so there is no need to build another waste treatment plant. The value products obtained through professional processing are sold separately, so the outflow logistics of the treatment plant are not studied, but only the simple reverse logistics process is studied here. Assumption 3. The whole recycling path of waste products can be roughly described as: recycling point-storage point-professional disposal point. The recycling operation period is one year, and the waste products produced during the recycling period will be recycled through this process. Assumption 4. This paper studies the product recycling treatment in the province, and the default transportation mode is road transportation. Therefore, the unit transportation cost between each facility is fixed, and the transportation cost is only related to the distance and transportation volume. Assumption 5. The fixed costs of building the same facilities should use uniform market prices, which do not vary from region to region. During the recycling period, no transportation activities are carried out between logistics facilities with the same function. Assumption 6. There is a quantitative limit on the amount of recycling accommodated by storage points and professional disposal points. However, the second-hand product market, component manufacturers and recycling material recyclers have no demand restriction on the recycled waste product flowing to them; that is, the amount of recycling of existing products has not reached market saturation.

Assumption 7.
Freight charges between different facilities in the same location are ignored by default. Assumption 8. Classified storage points have fixed unit storage costs, and professional disposal points have fixed operating costs, which are positively correlated with the amount of storage and processing, respectively.

•
Logistics facility I is the collection of recycling points, in the model, which is all the cities participating in product recycling, i ∈ I; P is the set of established detection classification points, p ∈ P; P min and P max are respectively the lower limit and upper limit of the number of detection classification points; D is the set of established collection processing points, d ∈ D; D min and D max are respectively the lower limit and upper limit of the number of collection processing points. s i is the amount of recycling at the recycling point i; S p is the amount of storage at the detection classification point P; S d is the amount of disposal at the recycling disposal point d.
l ip is the route distance between recycling point i and detection classification point p; l pd is the route distance between detection classification point p and recycling processing point d; M p and M d are respectively the volume limits of the detection classification point and the recycling processing point. Meanwhile, it is assumed that there is no volume limit when the product seller acts as the recycling point. X ip and X pd are respectively the transportation volume of waste products from the recycling point i to the detection classification point p and the transportation volume of waste products from the detection classification point p to the recycling and treatment point d.

•
Classification of waste products According to the availability of waste products, the products are divided into dismountable type c 1 , recyclable parts type c 2 , and recyclable parts type c 3 , and the proportion of the three types in the total amount of recycled products is l 1 , l 2 and l 3 . R 1 and R 2 are respectively the weight ratio of usable materials of dismountable products and usable parts and components of recyclable products. r 1 is the unit value of disassembled products after material extraction, and the unit is CNY/kg; r 2 is the market value of disassembled parts of recyclable parts, and the unit is CNY/kg; r 3 is the second-hand market value of products that can be renovated and sold for a second time, and the unit is CNY/unit. W is the average weight of the product.

•
Cost of recycling F i , F p and F d represent the initial construction costs of the recycling point, detection classification point and recycling treatment point, respectively, in which the construction cost of the recycling point is zero or ignored. y represents the transportation cost of waste products per unit distance per unit quantity, and the unit is CNY/kg·km. In the process of recycling, N p refers to the cost generated by the processing unit product at the detection classification point, and the unit is CNY/part. The recycling cost generated by the recycling point is determined by the type of recycling. The recycling cost for the disassembled product, the recyclable part and the product sold for the second time that can be repaired are respectively N 1 , N 2 and N 3 . The recycling price of waste products is q, and the unit is CNY/part. w y represents the carbon emission from waste product transportation per unit distance per unit quantity; w f represents the carbon emission from waste product processing per unit quantity at the detection classification center; w n represents carbon emission from the second refurbishment treatment in the recycling and treatment of waste products per unit quantity; w l represents carbon emission from the disassembly of recyclable parts in the recycling and treatment of waste products per unit quantity; w c represents carbon emission from the fine disassembly in the recycling and treatment of waste products per unit quantity; w o represents carbon emission from waste disposal of waste products per unit quantity.

•
Decision variable R p is the variable 0 or 1, indicating whether the detection classification point is established in city p. If it is 1, it means that this point is selected as the classification point; if it is 0, it means that no classification point has been established at this point. R d is a variable of 0 or 1, indicating whether a recycling processing point is established in city d. If it is 1, it means that this point is selected as the processing point; if it is 0, no processing point is established at this point. X ip is the quantity of waste products transported from recycling point i to detection classification point p, and X pd is the quantity of waste products transported from detection classification point p to recycling and processing point d.

Carbon Emission Cost Calculation
The purpose of reverse logistics is to save energy and reduce emissions. In the research boom of the low-carbon economy and green logistics, greenhouse gas emissions have always been the focus of attention [17][18][19], and the setting of logistics nodes is closely related to carbon emissions. How to integrate the idea of emission reduction into the system design and build a green logistics network plays a practical role. Some results have been achieved by taking carbon emissions as a function parameter and participating in the optimization process with a quantitative method [20][21][22][23][24]. Considering the carbon emission sources of a reverse logistics network, they are mainly divided into transport carbon emission and carbon emission processing [25].
The carbon emission in logistics transportation includes the indirect or direct carbon dioxide emission from various substances and vehicles [26]. According to the theory of energy consumption, vehicle carbon emissions are related to route slope, body weight and route length [27][28][29]. There are obvious uncertainties between people and the road environment, which are treated as deterministic factors in this paper. When the vehicle itself is used as a means of transport, it is purchased uniformly by the company, and the factors affecting vehicle carbon emissions such as engine and shape have also been fixed. Therefore, as long as carbon emissions are under fixed conditions, considering vehicle fuel consumption is directly linked to vehicle carbon emission quantification. Therefore, it can be summarized as: carbon emissions = fuel consumption * CO 2 emission coefficient [30]. Therefore, transport carbon emissions can be quantified by the carbon emission coefficient of recycled products per unit distance per unit weight [31]. To be specific, the transportation is divided into two stages: recycling center to sorting center and sorting center to processing center, and the carbon emissions are represented by distance and transportation volume.
Waste products also produce carbon emissions in the processing process, and for different processing processes, carbon emissions are different. The product treatment process can be divided into classification treatment at the detection classification center, secondary renovation treatment in recycling treatment, recyclable parts dismantling in recycling treatment, fine dismantling in recycling treatment and waste treatment. According to the different process carbon emission standards [32,33], the carbon emission of product recycling process can be obtained.

Establishment of Dual Objective Functions
To minimize carbon emissions and increase recycling benefits, a dual-objective facility location optimization model constrained by multiple resources, e.g., recycling volume and recycling facilities, was established [34,35].
Minimum carbon emission is where E 1 represents the total carbon emission in the whole recycling process, T 1 represents the carbon emission in the transportation process, and T 2 represents the carbon emission in the treatment process. For transport carbon emissions: For dealing with carbon emissions: Appl. Sci. 2021, 11, 6466 6 of 25 Maximum recycling benefits are where E 2 represents the recovery income net profit of waste products, c 1 represents recycling income, and c 2 represents the transportation cost of waste products, which is related to the transportation volume and transportation distance. c 3 represents product recovery cost and processing cost, where the processing cost is related to the processing stage and has different processing cost standards. c 4 represents the construction cost of the detection classification center and the recycling treatment center. For recycling income: For transportation costs: For processing costs: For the construction cost: To sum up, the model objective function can be expressed as: • Constraints Recycling amount of the recycling center is equal to the storage amount of the classification detection center and is also equal to the processing amount of the recycling treatment center: The transport volume from the recycling center to the classification detection center is equal to the recycling volume of the recycling center; the transport volume from the classification detection center to the recycling treatment center is equal to the storage amount of the classification detection center:

Gravity Algorithm
The gravity algorithm regards the optimal solution set as a space, and the single solution in the solution set is the scattered particles in this space [36]. These particles are attracted by each other and produce an aggregation effect. The cleverness of this algorithm lies in the correlation between the fitness of the solution and the gravity of the object [37]. Additionally, the particle with the lower gravity is attracted to the particle with the higher gravity. In the end, the space appears to be the result of the particles being clustered in an optimal position.
Assume that the total number of particles in the gravitational space is N, and the position of the ith particle is defined as: where x d i represents the position of particle i in the d dimension. Because of gravity, the particles move with varying accelerations, thus renewing their positions. The particle position update formula is as follows: where F d i represents the sum of the attractive forces on particle i in the d dimension; a d i (t) represents the d dimensional acceleration of particle i at time t; M ii (t) represents the inertia of particle i, and it is proportional to the mass of particles.
In the simplified model, the particle inertial mass can be expressed by the fitness value of the particle solution. The larger the fitness value is, the larger the inertia mass of the particle is, and the smaller the motion acceleration of the particle is, the more stable the particle position is. The specific formula of particle mass expressed by the solution fitness value is as follows: where f it i represents the fitness value of the particle solution; worst(t) represents the minimum fitness value of the solution set, representing the worst particle; best(t) represents the maximum fitness of the solution set, representing the optimal particle. The inertia value of the particle can be obtained from the particle mass: According to the law of universal gravitation, the attraction of a particle is proportional to its own gravity and the gravity of other particles. Therefore, the core formula of GSA algorithm-namely, the particle gravity formula-can be obtained: where F d ij represents the attraction of particle i attracted by particle j in the d dimension space; G(t) represents gravity constant, which is the correction constant used to calculate the gravity of particles. Its value decreases with time, indicating that particles tend to be stable; G 0 represents initial gravitational constant; T represents time constant, which refers to the maximum number of iterations in the algorithm; R ij (t) represents the distance of particle i and j at time t; α, ε represents the correction constant.

The Gravitational Particle Swarm Optimization Algorithm
Bringing in the idea of particle swarm optimization algorithm, by combining the velocity update formula and the universal gravitation algorithm, search speed and accuracy of the universal gravitation algorithm are improved [38]. Combining the advantages of both, Ma et al. proposed a novel hybrid optimization algorithm. It can well improve the algorithm performance [39]. Therefore, altering the core formula of the particle swarm optimization-particle velocity update formula, the velocity formula applied to universal gravitation algorithm is obtained.
To balance the two algorithms [40] and to make their respective advantages complementary to each other, the optimal value operator, p best , in particle swarm velocity formula is deleted. In addition, a balance between optimizing in the whole situation and searching for local optimization is sought. After referring to previous studies [41], the learning factors c 1 , c 2 in the velocity formula are improved. where η represents operator constant.
The specific steps of the algorithm are as follows: Step 1: Initialize the position and velocity of particles.
Step 2: Calculate the fitness values of the particle solution sets.
Step 3: Knowing the position of particles, the initial mass, initial force and initial acceleration of the particle can be obtained according to the Equations (18)- (24).
Step 4: After arriving at the optimal solution of the particle swarm, it will be considered as the motor direction of the particle.
Step 6: Judge whether or not the number of iterations or the result of the population reach the standard. If optimal solution sets were found, stop iterating. If not, go back to step 2 again.
Step 7: Stop cycling and output the optimal solution sets.

The Framework of the Proposed Multi-Objective Algorithm
Based on the above improvement of the universal gravitation algorithm, the hybrid gravitational particle swarm optimization algorithm is obtained. To apply this algorithm to the multi-objective model, there are five missions that need to be finished: individual (solution) coding, initial particle swarm formation, adaptive value calculation and Pareto solution screening, individual update and update of the dominant solution set.

•
Coding mode The coding form of the solution directly affects the efficiency of the algorithm. For this reason, this paper designs a coding scheme based on multi-constraints, which has lower decoding complexity and can easily perform mutation and crossover operations. According to the feature of questions, a solution is coded as an integer value vector, which contains two pieces of code. One section codes the location of the classification detection points. The other section codes the location of the recycling and processing center. As shown in Table 1, section X is the distribution of classification detection centers in urban agglomeration, and section Y is the distribution of recycling and processing centers in urban agglomeration. For example, X = {1, 1, 0, 1, 0, 0, 1}, the first element is 1, indicating that a classification detection center has been established in city 1, and the third element is 0, indicating that no classification detection center has been established in city 3, and so on.
The values of individuals in the initial particle swarm are randomly generated under multiple constraints. To obtain the initial individual, the particle is set as the length of individual 2N according to the particle coding mode, where N represents the number of urban agglomeration. For the constituent elements in an individual, each element has a state. If the element is a logistics facility setting point, its state is 1. If not, its state is 0. Usually, as the number of urban points increases, the number of site selection schemes for logistics facilities will also increase dramatically. Obviously, some of the site selection schemes generated initially are not feasible; that is, they cannot satisfy specific constraints. They need to be tweaked, and a common approach is to regenerate or rebuild. To make the algorithm quickly generate finite solutions and converge to feasible solutions, each individual in the initial population is initialized by Algorithm 1.

•
Adaptive value calculation and Pareto solution screening To evolve for the better, the particle swarm must be evaluated by comparing the fitness value of each individual (site selection scheme). This scheme contains two targets. A site selection scheme of the model can be decoded for a solution X. Decode X and use the target formula to calculate the result value of two targets. For the speed and accuracy of optimizing, the result value is transformed to some extent, and then the fitness value f 1 and f 2 related to the two resulting values are obtained.
After calculating fitness value, the next key mission is searching for the Pareto solution from the particle swarm. The Pareto optimal solution sets mean that no solution is equal to or superior to the Pareto optimal solution sets in every objective evaluation. The solutions are chosen by comparing every target to ensure that they are the Pareto optimal solution. Algorithm 2 describes how to obtain the Pareto solution by screening all the solutions.

Algorithm 2. Pareto Rank and Crowding Distance Calculation of Particle Swarm Individuals
Input: Pop (swarm) Output: P (record Pareto rank of the particle swarm); S (record particle crowding distance of the particle swarm) (1) Use the target formula to calculate the target value of particle swarm individuals. The transportation of waste products among logistics facilities is based on the principle of the nearest point. Only when the destination facility reaches its maximum capacity can they proceed to the next nearest point. The target value also needs to be transformed through fitness, and two targets value are transformed the fitness value with the smaller the better.

•
Individual update To evolve in an optimal direction for the whole particle swarm, the law of individual variation is crucial. The improved universal gravitation algorithm brings in the particle swarm velocity update formula. Figure 1 shows an example in which a particle individual uses the velocity update formula to finish its own evolution. Additionally, Algorithm 3 de-scribes how to obtain the updating solution of the individual particle by velocity updating.

Algorithm 3. Individual Update
Input: Pop, P, S Output: PopNew (New particle swarm) (1) i = 1 (2) While i < = NIND Do (3) Initialize the particle velocity V i = 0; (4) Calculate the particle accelerated velocity of individual i through the formula for calculating the acceleration of gravity. (5) The superior individuals in the dominant solution set were randomly selected as G Best . (6) Calculate the velocity V i of the particle i.   •

Update swarm
The non-dominated sorting algorithm is used to divide the sw tiple levels according to its Pareto rank. Meanwhile, an excellent store the top-level solutions found during the search. In every g dominated solutions in the current swarm are seen as candidate so excellent solution set. The Pareto rank of candidate solutions is hi distance is larger, so the solution is selected more easily and adde tion set.
• Improve Gravity multi-objective algorithm flow The multi-objective algorithm uses the number of maximum i mination criterion. That is, when the number of maximum iteratio put the Pareto solution and terminate the algorithm. Pareto optim lution is equal to or superior to the Pareto optimal solution sets in tion. The flow chart of this proposed algorithm is shown in Figure   1 1 1 0 0 1 1 0 0 1 1 1 0 0   1 1 1 0 0 1 0 0 1 1 1 1 0 0

•
Update swarm The non-dominated sorting algorithm is used to divide the swarm solution into multiple levels according to its Pareto rank. Meanwhile, an excellent solution set is used to store the top-level solutions found during the search. In every generation, all the non-dominated solutions in the current swarm are seen as candidate solutions of the updating excellent solution set. The Pareto rank of candidate solutions is higher, and the crowding distance is larger, so the solution is selected more easily and added to the excellent solution set.

•
Improve Gravity multi-objective algorithm flow The multi-objective algorithm uses the number of maximum iterations g max as the termination criterion. That is, when the number of maximum iterations g max is reached, output the Pareto solution and terminate the algorithm. Pareto optimality means that no solution is equal to or superior to the Pareto optimal solution sets in every objective evaluation. The flow chart of this proposed algorithm is shown in Figure 3.

Description of Basic Situation
A mobile phone recycling enterprise plans to set up a mobile phone recycling logi tics network in Jilin Province. The province is divided into nine mobile phone recyclin zones by cities, as shown in Figure 4. Through the analysis in the previous chapter, th enterprise takes the mobile phone sellers in each city as the recycling station, making eac city own a recycling station, where the recycling scope can cover most areas of the cit The enterprise plans to set up testing and classification centers and recycling and pro cessing centers from the nine recycling zones. In order to meet the demand for recyclin the number of testing and classification centers needs to be controlled between three an five, and the number of recycling and processing centers needs to be controlled betwee two and four.

Description of Basic Situation
A mobile phone recycling enterprise plans to set up a mobile phone recycling logistics network in Jilin Province. The province is divided into nine mobile phone recycling zones by cities, as shown in Figure 4. Through the analysis in the previous chapter, the enterprise takes the mobile phone sellers in each city as the recycling station, making each city own a recycling station, where the recycling scope can cover most areas of the city. The enterprise plans to set up testing and classification centers and recycling and processing centers from the nine recycling zones. In order to meet the demand for recycling, the number of testing and classification centers needs to be controlled between three and five, and the number of recycling and processing centers needs to be controlled between two and four. Technology. The frequency of mobile phone users replacing their mobile phones is affected by the upgrading of mobile phones and the contracts of mobile phone operators. Through market research and analysis, the mainstream mobile phone replacement strategy in the present mobile phone market is to replace every two years. Combined with data from the Ministry of Industry and Information Technology, this paper stipulates that each user owns 1.1 mobile phones on average.
From the update speed of mobile phones, nearly half of China's mobile phones are facing elimination every year. However, because the market of mobile phone recycling is not sound yet, the recycling rate of mobile phones has been kept at a low level. According to investigation, the recovery rate of mobile phones in China is about 1-2%, which is a huge gap compared with the 30% recovery rate in developed countries. The big gap between the number of discarded mobile phones and the amount of recycling shows that there is a lot of room for improvement in China's mobile phone recycling market. It is believed that with the development of a mobile phone recycling market, the recycling rate of used mobile phones will eventually reach a satisfactory level. Taking the reverse logistics facility established in this paper as an example, its perfect recycling process will greatly improve the recycling rate of mobile phones, so the preset recycling rate of mobile phones in Jilin Province would be 20%. Through calculation, the recycling situation of mobile phones in each city of Jilin Province is obtained, as shown in Table 2. At present, the recycling value of mobile phones is mainly divided into three parts. First, popular phones return to the market after cleaning and refurbishing. Second, defective popular phones are dismantled and sold as parts. Third, after professional processing and extraction of precious metals, the scrap machine can recover gold, silver, copper, palladium and other precious metals from the waste mobile phone [38]. By comparing popular phones in the market, it is assumed that the average weight of mobile phones is 125 g/unit, and the proportion of disassembled mobile phones, recyclable parts mobile phones, re-sold mobile phones that can be renovated and scrapped by each testing and classification center is 20%, 30%, 40% and 10%, respectively.
To better measure the recycling value of a mobile phone, the recycling of waste mobile phones is divided into three levels, which include a refurbishing machine, parts dismantling machine and material dismantling machine. The parts are divided into two parts: the main board and the casing. The processing cost and selling price list of the repairable machine and parts are obtained as shown in Table 3. Table 3. Part of disposal costs and prices of used mobile phones (unit: CNY/unit).

Mobile Phone Motherboard Mobile Phone Shell
Processing cost 60 15 5 Selling price 500 34 20 For the dismantling material, by referring to the data and comparing with the actual situation, mobile phone motherboards and mobile phone shells accounted for 80% of the weight and 20% of mobile phone, respectively. The material disassembly ratio of the two modules is listed as shown in Tables 4 and 5. In addition, the processing cost and selling price list of various materials can be obtained as shown in Table 6.  Recycling facilities are divided into three types: recycling station, detection classification station and professional processing plant. For the recycling station, we can cooperate with mobile phone sellers, so the fixed cost can be ignored. As for the detection and classification station and professional processing factories, China's mobile phone recycling is not centralized at present, and the corresponding facilities are rarely built. Based on the experience of the facility, the following assumptions are obtained.
(1) The initial capital for the construction of a classified storage station is CNY 140,000, the daily classified storage cost of used mobile phones is 14 CNY/unit, and the fuel consumption of used mobile phones per unit volume is 41,000 J/kg.
(2) The professional treatment plant needs special equipment and capabilities, so the initial construction fund is CNY 1.8 million. The unit disposal cost of the recycling processing station for a scrap treatment part is 7 CNY/unit. Waste fuel consumption per unit volume of processing is 89,000 J/kg.
(3) The cost per unit distance in the transportation process of used mobile phones is 35 CNY/km·ton. The carbon emission from the transport of used mobile phones per unit distance per unit quantity is 13 KG/km·kg.
Referring to the prices of mobile phone recyclers such as Huishubao, the recycling price of used mobile phones mostly fluctuates around CNY 200. In this paper, the purchase price is set at 200 CNY/unit. In addition, the transportation cost of mobile phone recycling is related to the transportation distance and the transportation weight as seen above. For this reason, the route distances between recycling areas in Jilin Province are shown in Table 7.  Changchun  0  122  116  168  264  443  263  112  350  Jilin  122  0  234  279  281  317  262  231  454  Siping  116  234  0  283  259  554  372  83  465  Sonyuan  168  279  283  0  431  578  429  278  191  Tonghua  264  281  259  431  0  479  60  182  614  Yanbian  North Korea  Autonomous prefecture   443  317  554  578  479  0  416  529  760   Baishan  263  262  372  429  60  416  0  254  611  Liaoyuan  112  231  83  278  182  529  254  0  459  Baicheng  350  454  465  191  614  760  611  459  0 Considering the carbon emissions in the transport and processing stages, the carbon dioxide emission coefficient of gasoline used in the disposal stage is 2.9, while the carbon dioxide emission coefficient of electric energy used in the processing stage is 0.9. The carbon emission standard of each type of waste mobile phone is obtained by using the carbon emission coefficient in the processing stage, as shown in Tables 8 and 9. Table 8. Carbon emission from partial treatment of mobile phones (kg/unit).

Whole Mobile Phone Mobile Phone Motherboard Mobile Phone Shell
Unit carbon emission 1.9 13 Table 9. Carbon emission from material treatment (kg/unit).

Results and Analysis
Given the reverse logistics network model of used mobile phones established in the above two chapters and aiming at the two objective functions of carbon emission minimization and revenue maximization, the improved universal gravitation algorithm was used as the optimization tool. The method was implemented in Matlab R2014 and ran on a PC with Intel(R)Core(TM)i5 CPU(2.53GHz/4.00g RAM) with Windows 7 operating system. In this paper, we set the parameters as follows: g max = 200, PS = 100 and limit = PS/2. Considering that the determination of these values is not the primary study task, the specific contents are not presented in this paper.

Optimization Results
Using the algorithm model, the two objectives of minimum carbon emission and maximum revenue were considered to be solved separately, and two optimization results were obtained. Then, the optimal site selection results were obtained by solving the optimization again for the two objectives.

•
Consider only the case where carbon emissions are minimal After the optimization model of the single object gravitation algorithm, the proposed algorithm was executed in g max = 200 generation to achieve full convergence. The minimum amount of carbon emission was 2,579,108,513 kg, and the evolutionary results are shown in Figure 5. The optimal location of carbon emission can be obtained in Figure 6.

Results and Analysis
Given the reverse logistics network model of used mobi above two chapters and aiming at the two objective functions zation and revenue maximization, the improved universal gra as the optimization tool. The method was implemented in M with Intel(R)Core(TM)i5 CPU(2.53GHz/4.00g RAM) with W In this paper, we set the parameters as follows: gmax = 200, PS sidering that the determination of these values is not the prim contents are not presented in this paper.

Optimization Results
Using the algorithm model, the two objectives of mini maximum revenue were considered to be solved separately, a were obtained. Then, the optimal site selection results were o mization again for the two objectives.

•
Consider only the case where carbon emissions are mini After the optimization model of the single object gravitat algorithm was executed in gmax = 200 generation to achieve mum amount of carbon emission was 2,579,108,513 kg, and shown in Figure 5. The optimal location of carbon emission c  . 2021, 11, x FOR PEER REVIEW Figure 6. Optimal location of carbon emission. Figure 6 represents the optimal location of carbon emissio represent the classification and testing centers, and the hollow cling and processing centers. In this figure, the X-axis indicat and the Y-axis indicates the latitude of the city. As shown in Fi is: Corresponding the optimal solution to the selected city, t center can be established in Changchun, Jilin, Siping, Tonghu cycling and processing center can be established in Cha Liaoyuan. Consequently, the carbon emissions of the recycling The carbon emissions of each type are shown in Table 10. It can the transport carbon emissions account for the largest part of sions, accounting for more than 95%, which reminds us tha recycling network, we should pay more attention to the opti duce the transport carbon emissions.  Figure 6. Optimal location of carbon emission. Figure 6 represents the optimal location of carbon emission, in which the solid points represent the classification and testing centers, and the hollow points represent the recycling and processing centers. In this figure, the X-axis indicates the longitude of the city, and the Y-axis indicates the latitude of the city. As shown in Figure 6, the optimal solution is: Corresponding the optimal solution to the selected city, the classification and testing center can be established in Changchun, Jilin, Siping, Tonghua and Baicheng, and the recycling and processing center can be established in Changchun, Jilin, Siping and Liaoyuan. Consequently, the carbon emissions of the recycling network can be minimized. The carbon emissions of each type are shown in Table 10. It can be seen from the table that the transport carbon emissions account for the largest part of the recycling carbon emissions, accounting for more than 95%, which reminds us that in the construction of the recycling network, we should pay more attention to the optimization of the route to reduce the transport carbon emissions. • Consider only the case of maximum benefit Through the optimization model of the single object gravitation algorithm, the proposed algorithm was executed in g max = 200 generation to achieve full convergence. The maximum value of the recovered income was CNY 206,683,179.8. Figure 7 shows the evolution result of the recovery revenue target. The optimal site selection of recovery income can be obtained in Figure 8.      Figure 8 represents the optimal site selection of recovery income, in which the solid points represent the classification and testing centers, and the hollow points represent the recycling and processing centers. In this figure, the X-axis indicates the longitude of the city, and the Y-axis indicates the latitude of the city. As shown in Figure 8, the optimal solution is: Corresponding the optimal solution to the selected city, the classification and testing center was established in Changchun, Jilin, Songyuan and Baishan, and the recycling and processing center was established in Jilin and Siping. At this time, the profit value of the recycling network could reach the maximum. The various benefit/cost values are shown in Table 11. It can be seen from the table that in the recovery cost, the recovery and processing cost accounts for the largest part, accounting for more than 90% of the total cost. In order to save cost and improve revenue, it is necessary to make a breakthrough in the recovery process. In addition, transportation cost is also a point that needs to be paid attention to. Reasonable planning of routes can save costs. • Consider dual targets of carbon emissions and recycling revenue It can be seen from the single objective analysis above that the carbon emission target is mainly affected by transportation, so the appropriate increase of classification detection stations and recycling and processing stations can effectively reduce carbon emissions, while the recycling income is mainly affected by the processing cost, so the recycling task can be completed with fewer logistics facilities to get more benefits. To sum up, in order to achieve the dual goals of maximizing benefits and minimizing carbon emissions, it is necessary to balance the number and location of logistics facilities. Figure 9 shows the convergence record of logistics site selection in the process of dual-objective evolution. It can be seen that after 100 times of evolution, the algorithm completed the balance of the dual-objective.  It can be seen from the single objective analysis above that the carbon emission targe is mainly affected by transportation, so the appropriate increase of classification detection stations and recycling and processing stations can effectively reduce carbon emissions while the recycling income is mainly affected by the processing cost, so the recycling task can be completed with fewer logistics facilities to get more benefits. To sum up, in orde to achieve the dual goals of maximizing benefits and minimizing carbon emissions, it i necessary to balance the number and location of logistics facilities. Figure 9 shows th convergence record of logistics site selection in the process of dual-objective evolution. I can be seen that after 100 times of evolution, the algorithm completed the balance of th dual-objective. The optimized Pareto solution set is partially shown in Table 12, and its binocula values are shown in Table 13, where different solutions represent the different selected The optimized Pareto solution set is partially shown in Table 12, and its binocular  values are shown in Table 13, where different solutions represent the different selected cities. Table 14 shows the changes in the selection of classification detection stations and recycling processing stations. Table 12. Two-objective partial Pareto solution.

Sensitivity Analysis
In order to test the sensitivity of the improved universal gravitation algorithm to the parameters of the model, the recovery rate and recovery price were taken as the analysis quantities. For the recovery rate, four cases of 10% reduction, 5% reduction, 5% increase and 10% increase were simulated on the basis of the original setting, which were compared with the original results, as shown in Table 15. The results show that when the recovery rate decreases, the site selection of recycling facilities basically remains unchanged. When the recovery rate of each recovery station increases by 5% and 10%, the possibility of alternative city 6 as the facility site is greatly increased. Through comparison, it can be seen that once the recovery rate changes, the site selection of recycling facilities will also change accordingly. However, the existing cities that have established logistics facilities have little change, and the main change is the construction of logistics facilities in remote areas.
To test the sensitivity of the improved gravitation algorithm to the parameters of the model, the facility capacity was taken as the analysis quantity. According to the facility capacity, four scenarios of 20% reduction, 40% reduction, 20% increase and 40% increase were simulated on the basis of the original setting, which were compared with the original results, as shown in Table 16. The results show that when the capacity of recycling facilities is reduced, the number of cities with recycling facilities will increase significantly, which brings about the change of transportation routes. When the capacity of each recycling facility increases by 20% and 40%, the number of recycling facilities will be reduced and fixed gradually, and transportation carbon emissions and transportation costs will also have corresponding changes. In conclusion, although the change of the capacity of recycling facilities will also have an impact on the site selection of recycling facilities, the recycling benefits will not have a large fluctuation.

Discussion
Based on the results of Section 4.2, the conclusions can be summarized by stating that the formulated reverse logistics facility location model and the proposed multi-objective algorithm are effective for solving the problem of reverse logistics location based on energy consumption.
To illustrate the advantages of this study, the contributions and significance can be summarized as follows. From the theoretical viewpoint, a double objective function optimization model of reverse logistics facility location was constructed in which the capacity of recycling center, classification center, processing center and the transportation volume between different recycling facilities were taken as constraints, and the two objective functions of low-carbon and maximum revenue were used to optimize the location scheme of logistics facilities. From the practical viewpoint, the aim of the multi-objective optimization of reverse logistics location is to put the findings into production and use in daily life. In addition, this study proposes the solution framework of the multi-objective problem of reverse logistics location, which can be used in different cities/areas with different features/constraints. This study presents an accurate and systematic method/tool to solve the multiobjective optimization problem of reverse logistics location and can assist researchers in better comprehending the multi-objective optimization theoretically, as well as assisting designers in developing a better security system in urban planning.

Conclusions
With the deepening of the research on resource reuse, product recycling has attracted the attention of enterprises and academic researchers. In the research on product recycling, the research on reverse logistics networks, especially the location of reverse logistics facilities, has always been the focus. In this paper, energy consumption was taken as a starting point, the carbon emissions and benefits in the reverse logistics network were taken as the functional objectives for modeling and optimization, and an improved gravity optimization algorithm was used to solve the reverse logistics facility location model, arriving at an excellent set of location schemes. Considering the influencing factors of network facilities construction, this paper evaluated the alternatives with multiple indexes, and made the final choice of facilities location. The main conclusions can be summarized as follows: (1) A double objective function optimization model of reverse logistics facility location was constructed. The three-level recycling facilities of recycling center, classification center and processing center were combined with the original logistics system of manufacturers, suppliers and consumers to form a closed-loop logistics system. In order to study the facility location problem of three-level reverse logistics, the capacity of recycling center, classification center, processing center and the transportation volume between different recycling facilities were taken as constraints, and the two objective functions of low-carbon and maximum revenue were used to optimize the location scheme of logistics facilities. The carbon emission target was divided into two parts: transportation carbon emission and treatment carbon emission. According to the different stages of transportation and treatment, the carbon emission in the process of product recovery was quantified, and the income was divided into four parts according to the different stages of recovery, and the total recovery income was calculated, respectively.
(2) The improved universal gravitation algorithm was used to solve the location model of reverse logistics facilities, and the Pareto set of location scheme was obtained. By combining the optimization characteristics of classical particle swarm optimization algorithm and universal gravitation algorithm, the two algorithms complemented each other, and a hybrid particle swarm optimization algorithm with better optimization performance was obtained. Taking the mobile phone recycling in Jilin Province as an example, it was proved that the hybrid algorithm has excellent convergence performance, meets the needs of solving the facility location model, and obtains the optimal solution set satisfying the double objectives. In addition, the accuracy and stability of the algorithm were verified by different recovery rate and facility capacity.
In future research, our work will focus on the following: (1) applying additional effort to real-time advanced methods for selecting Pareto solutions [42][43][44] and (2) after noting that the raw data have uncertain and imprecise features, integrating an advanced optimization method in the model [45,46].