# Omni-Channel Product Distribution Network Design by Using the Improved Particle Swarm Optimization Algorithm

Authors: Si Zhang; Heying Zhu; Xia Li; Yan Wang
Journal: Discrete Dynamics in Nature and Society
 (2019)

---

## Abstract

Distribution network under Omni-Channel integration contains many levels. There are one or more dealers at each level which forms a many-to-many distribution network. Consumers purchase a wide variety of products and their demands are uncertain, which constitutes a complex demand network and increases the complexity of the supply chain network. This paper focuses on the integrated optimization of supply chain distribution network and demand network and constructs the joint randomization planning model of location and routing. The goal is to minimize the total costs of the supply chain network under uncertain customer demands. Based on the traditional particle swarm optimization (PSO), this study introduces the collaborative idea to reduce the coding dimension, improves the boundary processing strategy, and adopts the mutation operator to expand the search space. A case study of distribution under Omni-Channel integration in a large enterprise was done. The validity of the model and the effectiveness of the proposed method were verified by numerical experiments.


---

## Body

## 1. Introduction

Over the last two decades, the product transportation of the retail industry was largely driven by the Internet Technology development. More and more retailers not only have the traditional offline stores but also start to have the online stores to satisfy the demands of customers from the Internet. As the emergence of new online channels, retailers begin to transform from traditionally storefront-based retailers to multi-channel retailers ([1–3]) and consider the operation management on two channels. This gave birth to the concept of Omni-Channel integration.Omni-Channel integration is the deep integration of online services, offline experience, and modern logistics in new retail era. It makes the supply chain more complex. The operation goal of Omni-Channel integration supply chain is to meet customer demands with the lowest costs. Saghiri et al. [4] proposed that it is a challenging task to analyze the distribution strategy of Omni-Channel integration. The supply chain usually contains many levels and each level has many sellers, which leads to the complexity of the supply chain network. Besides, the geographical locations of the sellers are spread all over the country. These factors have a significant impact on transportation costs. In addition, customer demands are often uncertain which will increase the costs of supply chain. Therefore the optimization for the supply chain with Omni-Channel is a challenge task and it is necessary to have some efficient approaches to help retailers reduce their costs. For example, Su’ning, one of the largest retailers with Omni-Channel supply chain in China, adopted a distribution method which is city-centered distributed to the county and rural areas step by step and constructed a three-level logistics distribution system. This system is “National Logistics Center—Regional Logistics Center—City Logistics Center”. As of June 2017, it has covered 297 prefecture-level cities and above, with about 2079 authorized chain stores and 833 authorized service outlets. It launched “100-City Half-Day” service to require products that are delivered to the customer within half day after the order is placed, which poses greater challenges to the distribution of logistics.Therefore, facing such a large supply chain network, building an efficient, high-quality, and low-cost distribution network is the key to improve competitiveness and sustainable development of enterprises. The main aim of this paper is to construct a highly efficient supply chain system to minimize costs and improve service quality.This paper is organized as follows. Section1 introduces the background of Omni-Channel distribution. The related literature is reviewed in Section 2 from the perspective of LRP, uncertain demand, and solving methods, respectively. In Section 3, a distribution network for Omni-Channel supply chain is proposed and a model for this problem is built. To solve this model, one algorithm based on the traditional PSO is developed in Section 4 by introducing the collaborative idea to reduce the coding dimension, improving the boundary processing strategy, and introducing the mutation operator to expand the search space. Sections 5 and 6 are the numerical results and sensitivity analysis based on one case study. Section 7 concludes the whole paper.


## 2. Literature Review

This paper mainly focuses on the location and routing problem (LRP) for Omni-Channel supply chain. Early LRP issues were limited to the location of a single facility. Maranzana [5] studied the problem of supplier location with the goal of minimizing transportation costs. Subsequently, many researchers conducted further research on the systematic and networking nature of transportation and location decisions. Govindan et al. [6] introduced a two-tier supply chain with a time window based on the supply chain network distribution of perishable food, which determined the number of distribution centers and facility locations. Abdulrazik et al. [7] optimized the tertiary supply chain under the background of producing oil palm in Malaysia. In addition, from the types of product distribution in the LRP, products transportation has expanded from single-product distribution to multi-product distribution. Geoffrion and Graves [8] developed a multi-product distribution model using Mixed-Integer Programming (MIP). Tiwari et al. [9] proposed an algorithm for solving multi-stage, multi-product supply chain network design. Viswanathan and Mathur [10] studied a distribution system consisting of a warehouse multiple products and multiple retailers. Arntzen [11] et al. used a mixed integer programming to establish a multi-product and multi-stage supply chain model for large multinational companies. Although research on the LRP changes from the initial single location planning to more consideration of the systemic and networking feature of the supply chain, the current researches on supply chain network are mainly focused on two-level or three-level system, and the distribution route planning is mostly one-to-many or many-to-one delivery. For the retail industry, it usually has many supply chain levels and many types of products in the distribution. Therefore, it is necessary to study the supply chain with more levels and more products. This paper is to design a five-level supply chain network for the LRP with multi-product distribution in retail industry.One another factor that cannot be ignored in Omni-Channel retail supply chain study is the uncertainties of customer demands. Different scholars used different methods to research uncertainty. Wang and He [12] introduced the uncertainty of demand and studied the problem of distribution route and warehouse location under the condition of uncertain demand. Zhen et al. [13] used the scenario method to solve the problem of uncertainty in the production process of automobiles. Tapiero and Soliman [14] applied optimal control theory to solve the problem of multi-product transportation, multi-regional production, and cross-time inventory planning under demand uncertainty. Baghalian et al. [15] used multi-objective optimization to model the location and routing problem in low-carbon supply chain and explored the impact of demand uncertainty. In general, customers have great differences in the demand for different products. It is necessary to form a personalized demand network. To describe customer demands more accurately, we set each customer's uncertain demand for different kinds of products to follow different normal distributions instead of letting all customers uniformly meet a normal distribution for all products uniformly.The LRP with multi-level, multi-product, and personalized demand is a very complex optimization problem. There are two main streams solving the large-scale optimization problem. The first stream is to use new evolutionary algorithms or add local search strategies (e.g. tabu search) to original algorithm. Zhang et al. [16] adopted tabu search algorithm to improve the accuracy and efficiency of cold chain product distribution. Javid and Azad [17] proposed a method of synchronous optimization of location and routing. Two-stage heuristic algorithm based on tabu search and simulated annealing was used to improve the solution space. Xia et al. [18] established the dual-objective programming model and designed an adaptive tabu search algorithm. The second stream of literature studies is to decompose the high-dimensional of large-scale complex problems into low-dimensional simple problems and optimize them separately. Hu et al. [19] studied the large-scale refrigerated truck scheduling problem and combined a variable neighborhood search with a particle swarm optimization to develop a two-stage decomposition algorithm. In 2009, Li and Yao [20] proposed a new co-evolutionary particle swarm optimization algorithm (CPSO), which added random grouping and adaptive weighting strategy to prove that it is effective in high-dimensional indivisible problems. The optimized 30-dimensional CPSO algorithm which has been proposed earlier can be extended to 1000-dimensional problems now. Compared with manufacture industry, distribution in retail industry usually have more supply chain levels, which results in the high dimensions of decision variables. We develop an improved PSO algorithm based on CPSO to efficiently solve the optimization problem.Therefore, in our paper, we aim to contribute to the current studies on this area from the following three aspects. Firstly, under the Omni-Channel distribution, the supply chain network level has been increased, and a many-to-many distribution network has been formed, which makes the supply chain network distribution more in line with the actual situation. Secondly, we allow each customer’s demands to follow different normal distributions for each product and form a demand network accordingly. The demand network can better predict customer demands. Thirdly, based on the traditional PSO, the algorithm introduces the collaborative idea to reduce the coding dimension, improves the boundary processing strategy, and introduces the mutation operator to expand the search space. This method improves the efficiency of the algorithm and provides a solution for solving high dimensional problems.


## 3. Model

### 3.1. Model Description

This model includes some common characteristics of distribution under Omni-Channel integration, so the conclusions of the model can be useful for many enterprises of supply-chain distribution under Omni-Channel integration. There are five levels in the entire supply chain. These levels include National Logistics Center A, Regional Logistics Center B, City Logistics Center C, Dealer D, and Customer E. Distribution network shall distribute m-type products to two types of customers. The first type of customers E1 is large customers such as companies; the demand of such customers is delivered directly by the national logistics center A. The second type of customers E2 is small supermarkets, like stores. The orders of such customers will be delivered through the regional logistics center B, the city logistics center C, and the dealer D. We can intuitively see the model in Figure1.Figure 1
Supply chain distribution network.According to the supply chain network, we built a mixed integer programming model to minimize the total supply chain costs. The total costs include three items: transportation costs of the products, fixed operating costs of logistics center, and carbon treatment costs which is to deal with carbon emission generated from the product transportation process.


### 3.2. The Index of the Model

a: The number of national logistics centers b: The number of regional logistics centers c: The number of city logistics centers d: The number of dealers e 1: The number of first type of customers e 2: The number of second type of customers m: The number of the category of products


### 3.3. Parameter Symbols

G a e 1 m: The unit transportation cost of product m delivered to customer e1 through national logistics center a G a b m: The unit transportation cost of product m delivered to regional logistics center b through national logistics center a G b c m: The unit transportation cost of product m delivered to city logistics center c through regional logistics center b G c d m: The unit transportation cost of product m delivered to dealer d through city logistics center c G d e 2 m: The unit transportation cost of product m delivered to customer e2 through dealer d d i s a e 1: The distance between national logistics center a and customer e1 d i s a b: The distance between national logistics center a and regional logistics center b d i s b c: The distance between regional logistics center b and city logistics center c d i s c d: The distance between city logistics center c and dealer d d i s d e 2: The distance between dealer d and customer e2 F a: The fixed operating costs of national logistics center located at point a F b: The fixed operating costs of regional logistics center located at point b F c: The fixed operating costs of city logistics center located at point c F d: The fixed operating costs of dealer located at point d O a m: The average inventory of product m in national logistics center located at point a R e 1 m: The demand of customer e1 for product m R e 2 m: The demand of customer e2 for product m T: Unit carbon treatment cost N: A large enough positive number


### 3.4. Decision Variables

Q a e 1 m: The quantity of product m from national logistics center a to customer e1 Q a b m: The quantity of product m from national logistics center a to regional logistics center b Q b c m: The quantity of product m from regional logistics center b to city logistics center c Q c d m: The quantity of product m from city logistics center c to dealer d Q d e 2 m: The quantity of product m from dealer d to customer e2 y a: Whether the national logistics center a operates, if it operates, then ya equals 1; otherwise the value is 0 y b: Whether the regional logistics center b operates, if it operates, then yb equals 1; otherwise the value is 0 y c: Whether the city logistics center c operates, if it operates, then yc equals 1; otherwise the value is 0 y d: Whether the dealer d operates, if it operates, then yd equals 1; otherwise the value is 0


### 3.5. Data Model

(1) m i n  f x (2) f x = ∑ A ∑ E 1 ∑ M G ae1 m Q a e 1 m d i s a e 1 + ∑ A ∑ B ∑ M G a b m Q a b m d i s a b + ∑ B ∑ C ∑ M G b c m Q b c m d i s b c + ∑ C ∑ D ∑ M G cd m Q c d m d i s c d + ∑ D ∑ E 2 ∑ M G de2 m Q de2 m d i s d e 2 + ∑ A F a y a + ∑ B F b y b + ∑ C F c y c + ∑ D F d y d + ∑ A ∑ E 1 ∑ M T Q a e 1 m dis ae1 + ∑ A ∑ B ∑ M T Q a b m d i s a b + ∑ B ∑ C ∑ M T Q b c m d i s b c + ∑ C ∑ D ∑ M T Q c d m d i s c d + ∑ D ∑ E 2 ∑ M T Q de2 m d i s d e 2 (3) ∑ A Q a e 1 m = R e 1 m  ∀ e 1 ∈ E 1 ;   ∀ m ∈ M (4) ∑ D Q d e 2 m = R e 2 m  ∀ e 2 ∈ E 2 ;   ∀ m ∈ M (5) ∑ A Q a b m ≥ ∑ C Q b c m  ∀ b ∈ B ;   ∀ m ∈ M (6) ∑ B Q b c m ≥ ∑ D Q c d m  ∀ c ∈ C ;   ∀ m ∈ M (7) ∑ C Q c d m ≥ ∑ E 2 Q d e 2 m  ∀ d ∈ D ;   ∀ m ∈ M (8) ∑ E 1 Q a e 1 m + ∑ B Q a b m ≤ O a m ∗ y a  ∀ a ∈ A ;   ∀ m ∈ M (9) ∑ A Q a b m - N × y b ≤ 0  ∀ b ∈ B ;   ∀ m ∈ M (10) ∑ B Q b c m - N × y c ≤ 0  ∀ c ∈ C ;   ∀ m ∈ M (11) ∑ C Q c d m - N × y d ≤ 0  ∀ d ∈ D ;   ∀ m ∈ M (12) y a , y b , y c , y d ∈ 0 , 1  ∀ a ∈ A ;   ∀ b ∈ B ;   ∀ c ∈ C ;   ∀ d ∈ D (13) Q a e 1 m , Q a b m , Q b c m , Q c d m , Q d e 2 m ≥ 0  ∀ a ∈ A ;   ∀ b ∈ B ;   ∀ c ∈ C ;   ∀ d ∈ D ;   ∀ e 1 ∈ E 1 ;   ∀ e 2 ∈ E 2The objective function (2) aims to minimize the total costs of the supply chain network. The constraint (3) requires that the demands of the customer e1 can only be directly delivered by the national logistics center a, and the constraint (4) states that the product can only be provided from the dealer d to the customer e2. Equations (5)-(7) specify that the volume of each logistics center is greater than or equal to the shipment at that point. Constraint (8) ensures that the inventory from the national logistics center a is greater than the total volume of shipped products, and the operation of the national logistics center a is represented by the 0-1 variable. Constraint (9) donates that if the regional logistics center b is not in operation, the quantity of product shipped to the regional logistics center b is zero. While constraint (10) and constraint (11) indicate that if city logistics center c and the dealer d do not operate, the volume of product shipped to the city logistics center c and the dealer d is zero. Constraint (12) and constraint (13) are the value constraints of the decision variable.


## 3.1. Model Description

This model includes some common characteristics of distribution under Omni-Channel integration, so the conclusions of the model can be useful for many enterprises of supply-chain distribution under Omni-Channel integration. There are five levels in the entire supply chain. These levels include National Logistics Center A, Regional Logistics Center B, City Logistics Center C, Dealer D, and Customer E. Distribution network shall distribute m-type products to two types of customers. The first type of customers E1 is large customers such as companies; the demand of such customers is delivered directly by the national logistics center A. The second type of customers E2 is small supermarkets, like stores. The orders of such customers will be delivered through the regional logistics center B, the city logistics center C, and the dealer D. We can intuitively see the model in Figure1.Figure 1
Supply chain distribution network.According to the supply chain network, we built a mixed integer programming model to minimize the total supply chain costs. The total costs include three items: transportation costs of the products, fixed operating costs of logistics center, and carbon treatment costs which is to deal with carbon emission generated from the product transportation process.


## 3.2. The Index of the Model

a: The number of national logistics centers b: The number of regional logistics centers c: The number of city logistics centers d: The number of dealers e 1: The number of first type of customers e 2: The number of second type of customers m: The number of the category of products


## 3.3. Parameter Symbols

G a e 1 m: The unit transportation cost of product m delivered to customer e1 through national logistics center a G a b m: The unit transportation cost of product m delivered to regional logistics center b through national logistics center a G b c m: The unit transportation cost of product m delivered to city logistics center c through regional logistics center b G c d m: The unit transportation cost of product m delivered to dealer d through city logistics center c G d e 2 m: The unit transportation cost of product m delivered to customer e2 through dealer d d i s a e 1: The distance between national logistics center a and customer e1 d i s a b: The distance between national logistics center a and regional logistics center b d i s b c: The distance between regional logistics center b and city logistics center c d i s c d: The distance between city logistics center c and dealer d d i s d e 2: The distance between dealer d and customer e2 F a: The fixed operating costs of national logistics center located at point a F b: The fixed operating costs of regional logistics center located at point b F c: The fixed operating costs of city logistics center located at point c F d: The fixed operating costs of dealer located at point d O a m: The average inventory of product m in national logistics center located at point a R e 1 m: The demand of customer e1 for product m R e 2 m: The demand of customer e2 for product m T: Unit carbon treatment cost N: A large enough positive number


## 3.4. Decision Variables

Q a e 1 m: The quantity of product m from national logistics center a to customer e1 Q a b m: The quantity of product m from national logistics center a to regional logistics center b Q b c m: The quantity of product m from regional logistics center b to city logistics center c Q c d m: The quantity of product m from city logistics center c to dealer d Q d e 2 m: The quantity of product m from dealer d to customer e2 y a: Whether the national logistics center a operates, if it operates, then ya equals 1; otherwise the value is 0 y b: Whether the regional logistics center b operates, if it operates, then yb equals 1; otherwise the value is 0 y c: Whether the city logistics center c operates, if it operates, then yc equals 1; otherwise the value is 0 y d: Whether the dealer d operates, if it operates, then yd equals 1; otherwise the value is 0


## 3.5. Data Model

(1) m i n  f x (2) f x = ∑ A ∑ E 1 ∑ M G ae1 m Q a e 1 m d i s a e 1 + ∑ A ∑ B ∑ M G a b m Q a b m d i s a b + ∑ B ∑ C ∑ M G b c m Q b c m d i s b c + ∑ C ∑ D ∑ M G cd m Q c d m d i s c d + ∑ D ∑ E 2 ∑ M G de2 m Q de2 m d i s d e 2 + ∑ A F a y a + ∑ B F b y b + ∑ C F c y c + ∑ D F d y d + ∑ A ∑ E 1 ∑ M T Q a e 1 m dis ae1 + ∑ A ∑ B ∑ M T Q a b m d i s a b + ∑ B ∑ C ∑ M T Q b c m d i s b c + ∑ C ∑ D ∑ M T Q c d m d i s c d + ∑ D ∑ E 2 ∑ M T Q de2 m d i s d e 2 (3) ∑ A Q a e 1 m = R e 1 m  ∀ e 1 ∈ E 1 ;   ∀ m ∈ M (4) ∑ D Q d e 2 m = R e 2 m  ∀ e 2 ∈ E 2 ;   ∀ m ∈ M (5) ∑ A Q a b m ≥ ∑ C Q b c m  ∀ b ∈ B ;   ∀ m ∈ M (6) ∑ B Q b c m ≥ ∑ D Q c d m  ∀ c ∈ C ;   ∀ m ∈ M (7) ∑ C Q c d m ≥ ∑ E 2 Q d e 2 m  ∀ d ∈ D ;   ∀ m ∈ M (8) ∑ E 1 Q a e 1 m + ∑ B Q a b m ≤ O a m ∗ y a  ∀ a ∈ A ;   ∀ m ∈ M (9) ∑ A Q a b m - N × y b ≤ 0  ∀ b ∈ B ;   ∀ m ∈ M (10) ∑ B Q b c m - N × y c ≤ 0  ∀ c ∈ C ;   ∀ m ∈ M (11) ∑ C Q c d m - N × y d ≤ 0  ∀ d ∈ D ;   ∀ m ∈ M (12) y a , y b , y c , y d ∈ 0 , 1  ∀ a ∈ A ;   ∀ b ∈ B ;   ∀ c ∈ C ;   ∀ d ∈ D (13) Q a e 1 m , Q a b m , Q b c m , Q c d m , Q d e 2 m ≥ 0  ∀ a ∈ A ;   ∀ b ∈ B ;   ∀ c ∈ C ;   ∀ d ∈ D ;   ∀ e 1 ∈ E 1 ;   ∀ e 2 ∈ E 2The objective function (2) aims to minimize the total costs of the supply chain network. The constraint (3) requires that the demands of the customer e1 can only be directly delivered by the national logistics center a, and the constraint (4) states that the product can only be provided from the dealer d to the customer e2. Equations (5)-(7) specify that the volume of each logistics center is greater than or equal to the shipment at that point. Constraint (8) ensures that the inventory from the national logistics center a is greater than the total volume of shipped products, and the operation of the national logistics center a is represented by the 0-1 variable. Constraint (9) donates that if the regional logistics center b is not in operation, the quantity of product shipped to the regional logistics center b is zero. While constraint (10) and constraint (11) indicate that if city logistics center c and the dealer d do not operate, the volume of product shipped to the city logistics center c and the dealer d is zero. Constraint (12) and constraint (13) are the value constraints of the decision variable.


## 4. Proposed Heuristic Method

The model of LRP presented in Section3 aims to minimize the total costs in product transportation process. In a large-scale transportation network, the code dimension is high, making it difficult to find an optimal solution in a limited time. To solve the problem, we resort to heuristic methods, such as particle swarm optimization (PSO). PSO is a stochastic, parallel optimization algorithm that does not require the properties of the optimized function to be divisible, steerable, and continuous. In addition, to further improve the accuracy and efficiency of the algorithm, the improved PSO is proposed.


### 4.1. Particle Swarm Optimization

PSO was proposed by Kennedy and Eberhart [21] and Lian [22] pointed out that this method was a stochastic optimization method based on swarm intelligence and PSO was inspired by the social behavior of bird flocking and their means of information exchange. Due to its easy implementation and fast convergence, PSO has been successfully applied to solve nonlinear optimization problems. Ai and Kachitvichyanukul [23], Goksal et al. [24], and Alinaghian et al. [25] designed an improved random topology particle swarm optimization algorithm (RT-PSO) for time-dependent vehicle routing problems. Hu et al. [19] developed a two-stage decomposition algorithm by combining variable neighborhood search algorithm and PSO to solve the large-scale refrigerated truck scheduling problem.In PSO, each particle represents a solution to the optimization problem. Particle initialization is generated in a random way, and each particle iteratively updates its position through personal best and global best. Each particle searches for the optimal fitness value in theD-dimensional feasible solution space. There are N particles in a group. In iteration t, the position of particle i is represented as a D-dimensional vector Xidt=xi1t,xi2t,…,xiDt,i=1,2,…,N. At the same time, the velocity of the particle i is also a D-dimensional vector, denoted as Vidt=(vi1t,vi2t,…,viDt),i=1,2,…,N. The optimal position searched by the particle i in iteration t is called the personal best, which is denoted as Pbest=pi1,pi2,…,piD,i=1,2,…,N. The optimal position searched by the group in iteration t is the global best, denoted as Gbest=pg1,pg2,…,pgD. After finding the two optimal values, each particle updates its velocity and position according to the following formula:(14)vidt+1=w∗vidt+c1∗r1∗Pbestidt-xidt+c2∗r2∗Gbestdt-xidt(15)xidt+1=xidt+vidt+1In the standard PSO, there are some shortcomings; for example,(1) when solving the high dimensional optimization problem, the increase of the dimension makes the search space expand into an exponential form, which is prone to premature convergence. That is, the particle group aggregates prematurely to the local optimal solution. (2) Particles may aggregate to local optimum too early in the process of gathering. If there is no escaping mechanism, the optimal solution may not be obtained. Therefore, we proposed improved strategy accordingly.


### 4.2. Strategies for Improving PSO

The improved PSO focuses on three aspects:(1) facing the “dimensional disaster” in large-scale problems, cooperative particle swarm optimization (CPSO) is proposed to transform large-scale complex problems into low-dimensional simple problems. (2) When particle flies out of the boundary, the boundary processing strategy is adopted to change the flight trajectory of the particles. (3) When the particles are aggregated, in order to prevent premature convergence to the local optimum, the variation factor is introduced to further explore the search space.


#### 4.2.1. CPSO for Dimensionality Reduction

As the dimension of the optimization problem increases, the performance of all algorithm drops dramatically. An effective solution is the cooperative evolution strategy. Bergh and Engelbrecht [26] combined cooperative evolution strategy with PSO and proposed cooperative particle swarm optimization (CPSO). The basic idea of CPSO is using K independent particle groups to search different dimension directions in the D-dimensional search space. Each particle swarm is updated independently, and no information is shared between the particle swarms. When calculating the fitness value, the position vectors of the optimal particles in the particle swarm are combined together to form a D-dimensional vector to calculate the fitness value.In this paper, to calculate the traffic volume between supply chain levels, the particle coding in standard PSO is 11-dimensionalX(A,B,C,D,E1,E2,m,A′,B′,C′,D′).The first 6 dimensions are the supply chain routing between neighboring level, the 7th dimension demonstrates different products, and the last 4 dimensions represent the node operation decision. According to CPSO, decision variable is classified by supply chain levels. Decision variable can be split into five types of 3-dimensional routing decision Q1A,E1,m,Q2A,B,m,Q3B,C,m,Q4C,D,m, Q5D,E2,m and four types of 1-dimensional node operation decision Y1A,Y2B,Y3C,Y4(D). The process is shown in Figure 2. Coding dimension is reduced and the algorithm speed is improved through CPSO.Figure 2
Collaborative dimension reduction.


#### 4.2.2. Improved Boundary Processing Strategy

In the search process, in order to improve the convergence efficiency and reduce the invalid search, the boundary conditions are usually set for the position and velocity to narrow the search range. In the traditional PSO, when a particle goes beyond the range, the general processing method is to set the particle on the boundary. That is,(16)IfXid<XminthenXid=Xmin(17)Else  ifXid>XmaxthenXid=XmaxThrough this method, after several iterations, a plurality of particles will aggregate toward the boundary in a plurality of dimensions and the flight path of the particles will tend to be the same. So the efficiency of the particle group is reduced. Thus, the following improvement strategies have been proposed:(18)IfXid<XminthenXid=Xmin+c3∗r3∗Xmin(19)Else  ifXid>XmaxthenXid=Xmax-c3∗r3∗Xmaxr 3 is a random number distributed between [0,1), and c3 is a constant on [0.01,0.5). The value of c3 will change according to the algorithm, the objective function, etc. This paper sets it to 0.05 to enhance the ability of the algorithm to jump out of the boundary region, so that the particles return to the feasible search space. After using the boundary strategy, on the one hand, it can ensure that the particles are located in the feasible search space. On the other hand, it can prevent the particles from accumulating too much to the boundary, resulting in local optimization on the boundary of the search space. It can maximize the search space and improve the quality of the solution. A comparison of the two strategies is shown in Figure 3.Figure 3
Improved boundary processing strategy.


#### 4.2.3. Mutation Operators to Expand Search Space

When PSO appears as local convergence or global convergence, the particles will have an aggregation phenomenon; that is, the particles have the same fitness. The fitness variance is the judgment of the convergence degree of the whole group. The larger the fitness variance is, the more likely the particle swarm is in the search state. When the fitness variance becomes smaller, the particle swarm is tending to converge; when the fitness variance is 0, PSO achieves local optimization or global optimization. The formula for calculating the fitness variance is as follows:(20)σ2=∑i=1mfi-favf2And,f=max1,max1≤i≤Nfi-fav, the function of f is to limit σ2.N is the total number of particles in the particle group.fi is the fitness value of particle i.fav=(1/N)∑i=1Nf is average fitness function value of each particle.If it is local optimum, the algorithm is premature. When the algorithm is premature, it will affect the accuracy of the algorithm. So,Gbestdt will be the personal best. The variation factor Pr is introduced, thereby changing the forward direction of the particle. So Gbestdt and Pbestidt will be updated:(21)IfPr>randdthen  update  particleThe probability of mutation is relatively small;Pr generally takes a number between [0, 1); in this algorithm, Pr=0.1;randd is a random number distributed between [0, 1). The specific approach is to mutate the particles when they are within 10% of the probability of mutation. Firstly, all the particles are sorted according to the fitness value; secondly, m particles with the best fitness value (m is usually half of the total number of particles) are obtained; then, the former Pr∗m particles are mutated as follows:(22)xidt+1=xidt∗1+0.5ηwhereη is a random variable obeying the standard normal distribution N(0,1). The purpose of using the mutation factor is to avoid the particle group falling into local optimum and better approaching the global optimum. When the algorithm appears premature, the mutation can change the direction of the entire particle swarm, and the particle swarm reaches a new field for search. By introducing the mutation operator, the algorithm can expand the search space when the search is stagnant. The improved PSO flow is shown in Figure 4. The pseudocode of improved PSO is presented in Box 1.Box 1:The improved PSO.
Initialization Set the parameters MaxItr,N,w,c1,c2,r1,r2,c3,r3,Pr; Initialize 9 independent particle swarms: Q a e 1 m a , e 1 , m , Q a b m a , b , m , Q b c m b , c , m , Q c d m c , d , m , Q d e 2 m d , e 2 , m y a a , y b b , y c c , y d ( d ) which are in the 11-dimensional solution space Xa,b,c,d,e1,e2,m,a′,b′,c′,d′; Define position (x) and velocity (v) for each particle; Assume individual best (Pbestidt) equal to the initialized particles; Select the best particle as global best (Gbestdt); Update Let 9 particle swarms be updated independently as follows: For t = 1 to MaxItr Update position (x) and velocity (v); v i d t + 1 = w ∗ v i d t + c 1 ∗ r 1 ∗ P b e s t i d t - x i d t + c 2 ∗ r 2 ∗ G b e s t d t - x i d t x i d t + 1 = x i d t + v i d t + 1 If  Xid<XminthenXid=Xmin+c3∗r3∗Xmin Else if  Xid>XmaxthenXid=Xmax-c3∗r3∗Xmax Respectively; End Evaluate the particles and calculate the fitness values; Sort the fitness values; If  Pr>randd x i d t + 1 = x i d t ∗ ( 1 + 0.5 η ) End Update personal best (Pbestidt) of 9 independent particle swarms; Update global best (Gbestdt) of 9 independent particle swarms; End Report Global best (Gbestdt) of 9 independent particle swarms; G b e s t d t found by splicing 9 independent particle swarms; Fitness values calculated by Gbestdt.Figure 4
Particle swarm optimization flowchart.


## 4.1. Particle Swarm Optimization

PSO was proposed by Kennedy and Eberhart [21] and Lian [22] pointed out that this method was a stochastic optimization method based on swarm intelligence and PSO was inspired by the social behavior of bird flocking and their means of information exchange. Due to its easy implementation and fast convergence, PSO has been successfully applied to solve nonlinear optimization problems. Ai and Kachitvichyanukul [23], Goksal et al. [24], and Alinaghian et al. [25] designed an improved random topology particle swarm optimization algorithm (RT-PSO) for time-dependent vehicle routing problems. Hu et al. [19] developed a two-stage decomposition algorithm by combining variable neighborhood search algorithm and PSO to solve the large-scale refrigerated truck scheduling problem.In PSO, each particle represents a solution to the optimization problem. Particle initialization is generated in a random way, and each particle iteratively updates its position through personal best and global best. Each particle searches for the optimal fitness value in theD-dimensional feasible solution space. There are N particles in a group. In iteration t, the position of particle i is represented as a D-dimensional vector Xidt=xi1t,xi2t,…,xiDt,i=1,2,…,N. At the same time, the velocity of the particle i is also a D-dimensional vector, denoted as Vidt=(vi1t,vi2t,…,viDt),i=1,2,…,N. The optimal position searched by the particle i in iteration t is called the personal best, which is denoted as Pbest=pi1,pi2,…,piD,i=1,2,…,N. The optimal position searched by the group in iteration t is the global best, denoted as Gbest=pg1,pg2,…,pgD. After finding the two optimal values, each particle updates its velocity and position according to the following formula:(14)vidt+1=w∗vidt+c1∗r1∗Pbestidt-xidt+c2∗r2∗Gbestdt-xidt(15)xidt+1=xidt+vidt+1In the standard PSO, there are some shortcomings; for example,(1) when solving the high dimensional optimization problem, the increase of the dimension makes the search space expand into an exponential form, which is prone to premature convergence. That is, the particle group aggregates prematurely to the local optimal solution. (2) Particles may aggregate to local optimum too early in the process of gathering. If there is no escaping mechanism, the optimal solution may not be obtained. Therefore, we proposed improved strategy accordingly.


## 4.2. Strategies for Improving PSO

The improved PSO focuses on three aspects:(1) facing the “dimensional disaster” in large-scale problems, cooperative particle swarm optimization (CPSO) is proposed to transform large-scale complex problems into low-dimensional simple problems. (2) When particle flies out of the boundary, the boundary processing strategy is adopted to change the flight trajectory of the particles. (3) When the particles are aggregated, in order to prevent premature convergence to the local optimum, the variation factor is introduced to further explore the search space.


### 4.2.1. CPSO for Dimensionality Reduction

As the dimension of the optimization problem increases, the performance of all algorithm drops dramatically. An effective solution is the cooperative evolution strategy. Bergh and Engelbrecht [26] combined cooperative evolution strategy with PSO and proposed cooperative particle swarm optimization (CPSO). The basic idea of CPSO is using K independent particle groups to search different dimension directions in the D-dimensional search space. Each particle swarm is updated independently, and no information is shared between the particle swarms. When calculating the fitness value, the position vectors of the optimal particles in the particle swarm are combined together to form a D-dimensional vector to calculate the fitness value.In this paper, to calculate the traffic volume between supply chain levels, the particle coding in standard PSO is 11-dimensionalX(A,B,C,D,E1,E2,m,A′,B′,C′,D′).The first 6 dimensions are the supply chain routing between neighboring level, the 7th dimension demonstrates different products, and the last 4 dimensions represent the node operation decision. According to CPSO, decision variable is classified by supply chain levels. Decision variable can be split into five types of 3-dimensional routing decision Q1A,E1,m,Q2A,B,m,Q3B,C,m,Q4C,D,m, Q5D,E2,m and four types of 1-dimensional node operation decision Y1A,Y2B,Y3C,Y4(D). The process is shown in Figure 2. Coding dimension is reduced and the algorithm speed is improved through CPSO.Figure 2
Collaborative dimension reduction.


### 4.2.2. Improved Boundary Processing Strategy

In the search process, in order to improve the convergence efficiency and reduce the invalid search, the boundary conditions are usually set for the position and velocity to narrow the search range. In the traditional PSO, when a particle goes beyond the range, the general processing method is to set the particle on the boundary. That is,(16)IfXid<XminthenXid=Xmin(17)Else  ifXid>XmaxthenXid=XmaxThrough this method, after several iterations, a plurality of particles will aggregate toward the boundary in a plurality of dimensions and the flight path of the particles will tend to be the same. So the efficiency of the particle group is reduced. Thus, the following improvement strategies have been proposed:(18)IfXid<XminthenXid=Xmin+c3∗r3∗Xmin(19)Else  ifXid>XmaxthenXid=Xmax-c3∗r3∗Xmaxr 3 is a random number distributed between [0,1), and c3 is a constant on [0.01,0.5). The value of c3 will change according to the algorithm, the objective function, etc. This paper sets it to 0.05 to enhance the ability of the algorithm to jump out of the boundary region, so that the particles return to the feasible search space. After using the boundary strategy, on the one hand, it can ensure that the particles are located in the feasible search space. On the other hand, it can prevent the particles from accumulating too much to the boundary, resulting in local optimization on the boundary of the search space. It can maximize the search space and improve the quality of the solution. A comparison of the two strategies is shown in Figure 3.Figure 3
Improved boundary processing strategy.


### 4.2.3. Mutation Operators to Expand Search Space

When PSO appears as local convergence or global convergence, the particles will have an aggregation phenomenon; that is, the particles have the same fitness. The fitness variance is the judgment of the convergence degree of the whole group. The larger the fitness variance is, the more likely the particle swarm is in the search state. When the fitness variance becomes smaller, the particle swarm is tending to converge; when the fitness variance is 0, PSO achieves local optimization or global optimization. The formula for calculating the fitness variance is as follows:(20)σ2=∑i=1mfi-favf2And,f=max1,max1≤i≤Nfi-fav, the function of f is to limit σ2.N is the total number of particles in the particle group.fi is the fitness value of particle i.fav=(1/N)∑i=1Nf is average fitness function value of each particle.If it is local optimum, the algorithm is premature. When the algorithm is premature, it will affect the accuracy of the algorithm. So,Gbestdt will be the personal best. The variation factor Pr is introduced, thereby changing the forward direction of the particle. So Gbestdt and Pbestidt will be updated:(21)IfPr>randdthen  update  particleThe probability of mutation is relatively small;Pr generally takes a number between [0, 1); in this algorithm, Pr=0.1;randd is a random number distributed between [0, 1). The specific approach is to mutate the particles when they are within 10% of the probability of mutation. Firstly, all the particles are sorted according to the fitness value; secondly, m particles with the best fitness value (m is usually half of the total number of particles) are obtained; then, the former Pr∗m particles are mutated as follows:(22)xidt+1=xidt∗1+0.5ηwhereη is a random variable obeying the standard normal distribution N(0,1). The purpose of using the mutation factor is to avoid the particle group falling into local optimum and better approaching the global optimum. When the algorithm appears premature, the mutation can change the direction of the entire particle swarm, and the particle swarm reaches a new field for search. By introducing the mutation operator, the algorithm can expand the search space when the search is stagnant. The improved PSO flow is shown in Figure 4. The pseudocode of improved PSO is presented in Box 1.Box 1:The improved PSO.
Initialization Set the parameters MaxItr,N,w,c1,c2,r1,r2,c3,r3,Pr; Initialize 9 independent particle swarms: Q a e 1 m a , e 1 , m , Q a b m a , b , m , Q b c m b , c , m , Q c d m c , d , m , Q d e 2 m d , e 2 , m y a a , y b b , y c c , y d ( d ) which are in the 11-dimensional solution space Xa,b,c,d,e1,e2,m,a′,b′,c′,d′; Define position (x) and velocity (v) for each particle; Assume individual best (Pbestidt) equal to the initialized particles; Select the best particle as global best (Gbestdt); Update Let 9 particle swarms be updated independently as follows: For t = 1 to MaxItr Update position (x) and velocity (v); v i d t + 1 = w ∗ v i d t + c 1 ∗ r 1 ∗ P b e s t i d t - x i d t + c 2 ∗ r 2 ∗ G b e s t d t - x i d t x i d t + 1 = x i d t + v i d t + 1 If  Xid<XminthenXid=Xmin+c3∗r3∗Xmin Else if  Xid>XmaxthenXid=Xmax-c3∗r3∗Xmax Respectively; End Evaluate the particles and calculate the fitness values; Sort the fitness values; If  Pr>randd x i d t + 1 = x i d t ∗ ( 1 + 0.5 η ) End Update personal best (Pbestidt) of 9 independent particle swarms; Update global best (Gbestdt) of 9 independent particle swarms; End Report Global best (Gbestdt) of 9 independent particle swarms; G b e s t d t found by splicing 9 independent particle swarms; Fitness values calculated by Gbestdt.Figure 4
Particle swarm optimization flowchart.


## 4.2.1. CPSO for Dimensionality Reduction

As the dimension of the optimization problem increases, the performance of all algorithm drops dramatically. An effective solution is the cooperative evolution strategy. Bergh and Engelbrecht [26] combined cooperative evolution strategy with PSO and proposed cooperative particle swarm optimization (CPSO). The basic idea of CPSO is using K independent particle groups to search different dimension directions in the D-dimensional search space. Each particle swarm is updated independently, and no information is shared between the particle swarms. When calculating the fitness value, the position vectors of the optimal particles in the particle swarm are combined together to form a D-dimensional vector to calculate the fitness value.In this paper, to calculate the traffic volume between supply chain levels, the particle coding in standard PSO is 11-dimensionalX(A,B,C,D,E1,E2,m,A′,B′,C′,D′).The first 6 dimensions are the supply chain routing between neighboring level, the 7th dimension demonstrates different products, and the last 4 dimensions represent the node operation decision. According to CPSO, decision variable is classified by supply chain levels. Decision variable can be split into five types of 3-dimensional routing decision Q1A,E1,m,Q2A,B,m,Q3B,C,m,Q4C,D,m, Q5D,E2,m and four types of 1-dimensional node operation decision Y1A,Y2B,Y3C,Y4(D). The process is shown in Figure 2. Coding dimension is reduced and the algorithm speed is improved through CPSO.Figure 2
Collaborative dimension reduction.


## 4.2.2. Improved Boundary Processing Strategy

In the search process, in order to improve the convergence efficiency and reduce the invalid search, the boundary conditions are usually set for the position and velocity to narrow the search range. In the traditional PSO, when a particle goes beyond the range, the general processing method is to set the particle on the boundary. That is,(16)IfXid<XminthenXid=Xmin(17)Else  ifXid>XmaxthenXid=XmaxThrough this method, after several iterations, a plurality of particles will aggregate toward the boundary in a plurality of dimensions and the flight path of the particles will tend to be the same. So the efficiency of the particle group is reduced. Thus, the following improvement strategies have been proposed:(18)IfXid<XminthenXid=Xmin+c3∗r3∗Xmin(19)Else  ifXid>XmaxthenXid=Xmax-c3∗r3∗Xmaxr 3 is a random number distributed between [0,1), and c3 is a constant on [0.01,0.5). The value of c3 will change according to the algorithm, the objective function, etc. This paper sets it to 0.05 to enhance the ability of the algorithm to jump out of the boundary region, so that the particles return to the feasible search space. After using the boundary strategy, on the one hand, it can ensure that the particles are located in the feasible search space. On the other hand, it can prevent the particles from accumulating too much to the boundary, resulting in local optimization on the boundary of the search space. It can maximize the search space and improve the quality of the solution. A comparison of the two strategies is shown in Figure 3.Figure 3
Improved boundary processing strategy.


## 4.2.3. Mutation Operators to Expand Search Space

When PSO appears as local convergence or global convergence, the particles will have an aggregation phenomenon; that is, the particles have the same fitness. The fitness variance is the judgment of the convergence degree of the whole group. The larger the fitness variance is, the more likely the particle swarm is in the search state. When the fitness variance becomes smaller, the particle swarm is tending to converge; when the fitness variance is 0, PSO achieves local optimization or global optimization. The formula for calculating the fitness variance is as follows:(20)σ2=∑i=1mfi-favf2And,f=max1,max1≤i≤Nfi-fav, the function of f is to limit σ2.N is the total number of particles in the particle group.fi is the fitness value of particle i.fav=(1/N)∑i=1Nf is average fitness function value of each particle.If it is local optimum, the algorithm is premature. When the algorithm is premature, it will affect the accuracy of the algorithm. So,Gbestdt will be the personal best. The variation factor Pr is introduced, thereby changing the forward direction of the particle. So Gbestdt and Pbestidt will be updated:(21)IfPr>randdthen  update  particleThe probability of mutation is relatively small;Pr generally takes a number between [0, 1); in this algorithm, Pr=0.1;randd is a random number distributed between [0, 1). The specific approach is to mutate the particles when they are within 10% of the probability of mutation. Firstly, all the particles are sorted according to the fitness value; secondly, m particles with the best fitness value (m is usually half of the total number of particles) are obtained; then, the former Pr∗m particles are mutated as follows:(22)xidt+1=xidt∗1+0.5ηwhereη is a random variable obeying the standard normal distribution N(0,1). The purpose of using the mutation factor is to avoid the particle group falling into local optimum and better approaching the global optimum. When the algorithm appears premature, the mutation can change the direction of the entire particle swarm, and the particle swarm reaches a new field for search. By introducing the mutation operator, the algorithm can expand the search space when the search is stagnant. The improved PSO flow is shown in Figure 4. The pseudocode of improved PSO is presented in Box 1.Box 1:The improved PSO.
Initialization Set the parameters MaxItr,N,w,c1,c2,r1,r2,c3,r3,Pr; Initialize 9 independent particle swarms: Q a e 1 m a , e 1 , m , Q a b m a , b , m , Q b c m b , c , m , Q c d m c , d , m , Q d e 2 m d , e 2 , m y a a , y b b , y c c , y d ( d ) which are in the 11-dimensional solution space Xa,b,c,d,e1,e2,m,a′,b′,c′,d′; Define position (x) and velocity (v) for each particle; Assume individual best (Pbestidt) equal to the initialized particles; Select the best particle as global best (Gbestdt); Update Let 9 particle swarms be updated independently as follows: For t = 1 to MaxItr Update position (x) and velocity (v); v i d t + 1 = w ∗ v i d t + c 1 ∗ r 1 ∗ P b e s t i d t - x i d t + c 2 ∗ r 2 ∗ G b e s t d t - x i d t x i d t + 1 = x i d t + v i d t + 1 If  Xid<XminthenXid=Xmin+c3∗r3∗Xmin Else if  Xid>XmaxthenXid=Xmax-c3∗r3∗Xmax Respectively; End Evaluate the particles and calculate the fitness values; Sort the fitness values; If  Pr>randd x i d t + 1 = x i d t ∗ ( 1 + 0.5 η ) End Update personal best (Pbestidt) of 9 independent particle swarms; Update global best (Gbestdt) of 9 independent particle swarms; End Report Global best (Gbestdt) of 9 independent particle swarms; G b e s t d t found by splicing 9 independent particle swarms; Fitness values calculated by Gbestdt.Figure 4
Particle swarm optimization flowchart.


## 5. Numerical Experiments

In this section, some numerical experiments are performed to verify the effectiveness of the proposed model of multi-product distribution and location problem based on uncertain demand networks. The programming language is C#, and the experiments are completed on a PC (Intel Core i7, 2.6 GHz; Memory, 8G).


### 5.1. Parameter Setting

According to the background of the distribution network under Omni-Channel integration, the parameters of this experiment are set. The location and routing process is the same as the model description in Section3. To demonstrate the broad applicability of the model, some of the parameters of the experiment were generated based on a random probability distribution shown in Table 1.Table 1
Parameters setting for the experiments.
Parameter Distribution type Distributions Unit G a e 1 m   G a b m   G b c m   G c d m   G d e 2 m Uniform distribution U ( 190,250 ) RMB/ton/km F a   F b   F c   F d Uniform distribution U ( 150000,350000 ) RMB/year O a m Uniform distribution U ( 150000,200000 ) ton R Normal distribution N ( μ j p , σ 2 ) ton /yearThe unit transportation costG is uniformly generated according to different road conditions. The operating costs F of dealers for different levels are uniformly generated according to the land price and scale. The capacity O of the distribution center is uniformly distributed according to the product category and scale. This article deals with the uncertain customer demand as the following two steps.Step 1.
Each customer has different demands for different products. A set of scenarios with different demands are set according to the characteristics of the normal distribution function.Step 2.
Let different customer demands for different products meet different normal distributions. This creates a network of demand of each customer for each product, rather than allowing all customers to meet only one normal distribution demand for each product, which could predict customer demands more accurately.


### 5.2. Experimental Results

This paper studies multi-product distribution under the uncertainty of customer demand, so the experiment is divided into four groups according to the changes of customer demand. Each group of experiments is divided into five different situations according to the number of suppliers and customers. The results in Table2 show 5 experimental scales with different changes of demand. The objective function value and running time are calculated by the improved PSO and CPLEX solver. “1-1-2-4-10 (2+8)” represents one national logistics center, one regional logistics center, two city logistics centers, four dealers, and ten customers (two of them are first-type customers and eight of them are second-type customers).Table 2
Comparison of different kinds of algorithms.
CPLEX Improved PSO Gap Scale of the study Obj. CPU(s) Obj. CPU(s) [(4)-(2)]/(2) ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 ) Change in demand =0 1-1-2-4-10 (2+8) 155118311.9 9.95 155118311.9 13.31 0.00% 2-2-4-10-35 (4+31) 4850885922 26.8 4861980932 29.03 0.23% 2-2-6-14-55 (5+50) 5181467988 45 5213955899 42 0.63% 2-4-12-22-85 (7+78) 5318725918 304.1 5363403216 77.06 0.84% 3-5-18-70-205 (10+195) 5619609011 >3600 5672995297 503 0.95%  Change in demand =0.01 1-1-2-4-10 (2+8) 155478693.3 9.95 155478693.3 11.71 0.00% 2-2-4-10-35 (4+31) 4905668616 23.45 4913027119 28.57 0.15% 2-2-6-14-55 (5+50) 5193458934 49 5224634879 44.54 0.60% 2-4-12-22-85 (7+78) 5304694312 325 5350314683 83.86 0.86% 3-5-18-70-205 (10+195) 5608465011 >3600 5666793047 524 1.04%  Change in demand =0.05 1-1-2-4-10 (2+8) 155446843.3 9.9 155446843.3 8.63 0.00% 2-2-4-10-35 (4+31) 4902814227 23.75 4913110137 24.77 0.21% 2-2-6-14-55 (5+50) 5220507224 45.73 5231902083 38.29 0.22% 2-4-12-22-85 (7+78) 5342902440 341.6 5392591433 76.8 0.93% 3-5-18-70-205 (10+195) 5621433450 >3600 5685517791 537 1.14%  Change in demand =0.1 1-1-2-4-10 (2+8) 155867650.3 9.95 155867650.3 12.71 0.00% 2-2-4-10-35 (4+31) 4928341567 29.9 4941648089 27.94 0.27% 2-2-6-14-55 (5+50) 5210996332 57.91 5233955899 48.36 0.44% 2-4-12-22-85 (7+78) 5339691338 364.7 5389350467 74.74 0.93% 3-5-18-70-205 (10+195) 5620433450 >3600 5691250911 556 1.26%  Average value      0.53%In small-scale cases, the running time of CPLEX is less than the improved PSO. As the scale increases, the running time of CPLEX increases significantly. For example, in the first three experiments, improved PSO and CPLEX have similar effects. When the scale increased to 85 customers, the running time of CPLEX increased quickly, while the running time of improved PSO is less increased. In general, the results of the CPLEX solver are considered to be the optimal solution. Therefore, the CPLEX solver solves the model faster than the improved PSO method on a small scale. However, when the network exceeds a certain scale, the improved PSO gets results faster than the CPLEX solver. In addition, the average deviation between the improved PSO and the optimal result is about 0.53%, which indicates that the result is accurate and the proposed algorithm is suitable for solving the model.


## 5.1. Parameter Setting

According to the background of the distribution network under Omni-Channel integration, the parameters of this experiment are set. The location and routing process is the same as the model description in Section3. To demonstrate the broad applicability of the model, some of the parameters of the experiment were generated based on a random probability distribution shown in Table 1.Table 1
Parameters setting for the experiments.
Parameter Distribution type Distributions Unit G a e 1 m   G a b m   G b c m   G c d m   G d e 2 m Uniform distribution U ( 190,250 ) RMB/ton/km F a   F b   F c   F d Uniform distribution U ( 150000,350000 ) RMB/year O a m Uniform distribution U ( 150000,200000 ) ton R Normal distribution N ( μ j p , σ 2 ) ton /yearThe unit transportation costG is uniformly generated according to different road conditions. The operating costs F of dealers for different levels are uniformly generated according to the land price and scale. The capacity O of the distribution center is uniformly distributed according to the product category and scale. This article deals with the uncertain customer demand as the following two steps.Step 1.
Each customer has different demands for different products. A set of scenarios with different demands are set according to the characteristics of the normal distribution function.Step 2.
Let different customer demands for different products meet different normal distributions. This creates a network of demand of each customer for each product, rather than allowing all customers to meet only one normal distribution demand for each product, which could predict customer demands more accurately.


## 5.2. Experimental Results

This paper studies multi-product distribution under the uncertainty of customer demand, so the experiment is divided into four groups according to the changes of customer demand. Each group of experiments is divided into five different situations according to the number of suppliers and customers. The results in Table2 show 5 experimental scales with different changes of demand. The objective function value and running time are calculated by the improved PSO and CPLEX solver. “1-1-2-4-10 (2+8)” represents one national logistics center, one regional logistics center, two city logistics centers, four dealers, and ten customers (two of them are first-type customers and eight of them are second-type customers).Table 2
Comparison of different kinds of algorithms.
CPLEX Improved PSO Gap Scale of the study Obj. CPU(s) Obj. CPU(s) [(4)-(2)]/(2) ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) ( 6 ) Change in demand =0 1-1-2-4-10 (2+8) 155118311.9 9.95 155118311.9 13.31 0.00% 2-2-4-10-35 (4+31) 4850885922 26.8 4861980932 29.03 0.23% 2-2-6-14-55 (5+50) 5181467988 45 5213955899 42 0.63% 2-4-12-22-85 (7+78) 5318725918 304.1 5363403216 77.06 0.84% 3-5-18-70-205 (10+195) 5619609011 >3600 5672995297 503 0.95%  Change in demand =0.01 1-1-2-4-10 (2+8) 155478693.3 9.95 155478693.3 11.71 0.00% 2-2-4-10-35 (4+31) 4905668616 23.45 4913027119 28.57 0.15% 2-2-6-14-55 (5+50) 5193458934 49 5224634879 44.54 0.60% 2-4-12-22-85 (7+78) 5304694312 325 5350314683 83.86 0.86% 3-5-18-70-205 (10+195) 5608465011 >3600 5666793047 524 1.04%  Change in demand =0.05 1-1-2-4-10 (2+8) 155446843.3 9.9 155446843.3 8.63 0.00% 2-2-4-10-35 (4+31) 4902814227 23.75 4913110137 24.77 0.21% 2-2-6-14-55 (5+50) 5220507224 45.73 5231902083 38.29 0.22% 2-4-12-22-85 (7+78) 5342902440 341.6 5392591433 76.8 0.93% 3-5-18-70-205 (10+195) 5621433450 >3600 5685517791 537 1.14%  Change in demand =0.1 1-1-2-4-10 (2+8) 155867650.3 9.95 155867650.3 12.71 0.00% 2-2-4-10-35 (4+31) 4928341567 29.9 4941648089 27.94 0.27% 2-2-6-14-55 (5+50) 5210996332 57.91 5233955899 48.36 0.44% 2-4-12-22-85 (7+78) 5339691338 364.7 5389350467 74.74 0.93% 3-5-18-70-205 (10+195) 5620433450 >3600 5691250911 556 1.26%  Average value      0.53%In small-scale cases, the running time of CPLEX is less than the improved PSO. As the scale increases, the running time of CPLEX increases significantly. For example, in the first three experiments, improved PSO and CPLEX have similar effects. When the scale increased to 85 customers, the running time of CPLEX increased quickly, while the running time of improved PSO is less increased. In general, the results of the CPLEX solver are considered to be the optimal solution. Therefore, the CPLEX solver solves the model faster than the improved PSO method on a small scale. However, when the network exceeds a certain scale, the improved PSO gets results faster than the CPLEX solver. In addition, the average deviation between the improved PSO and the optimal result is about 0.53%, which indicates that the result is accurate and the proposed algorithm is suitable for solving the model.


## 6. Case Analysis

Based on the background of the distribution of three types of products in Chengdu, China, the scale is “2-4-12-22-85(7+78)”. M Company is a typical retail enterprise and products transportation is under Omni-Channel integration. Products include electrical appliances, household products, and daily chemical products. The average annual demand of three products for 85 customers is shown in Table3. Their location is shown in Figure 5.Table 3
Annual average demand of 3 products for 85 customers.
No. D1 D2 D3 No. D1 D2 D3 No. D1 D2 D3 No. D1 D2 D3 1 99 9 12 23 41 27 3 45 24 4 5 67 29 5 1 2 21 7 1 24 44 1 1 46 14 4 4 68 10 2 2 3 68 18 3 25 6 5 2 47 4 2 2 69 33 8 1 4 4 2 2 26 2 0 1 48 15 12 2 70 56 9 2 5 55 2 1 27 104 23 1 49 19 7 4 71 4 1 3 6 16 6 3 28 12 4 2 50 2 1 0 72 43 23 2 7 14 5 1 29 1 1 1 51 9 7 0 73 16 5 2 8 3 1 0 30 1 1 1 52 4 1 1 74 9 1 1 9 17 3 1 31 5 2 1 53 16 7 7 75 24 2 3 10 1 0 0 32 8 1 1 54 14 8 1 76 2 1 1 11 4 1 7 33 48 9 7 55 22 7 2 77 7 3 2 12 18 1 6 34 23 15 1 56 4 3 1 78 23 5 2 13 67 14 7 35 65 8 4 57 6 2 2 79 163 19 1 14 3 1 2 36 31 15 1 58 1 1 2 80 429 4 1 15 9 5 1 37 31 3 4 59 4 3 0 81 417 54 1 16 7 3 1 38 9 1 1 60 13 7 3 82 501 25 1 17 43 3 15 39 16 2 1 61 13 3 1 83 286 22 4 18 1 0 7 40 45 11 2 62 12 2 2 84 455 25 2 19 29 2 3 41 86 8 1 63 45 8 3 85 393 45 2 20 28 8 4 42 39 5 4 64 39 14 3 21 38 6 1 43 16 3 1 65 23 6 1 22 5 2 3 44 6 2 1 66 67 7 1 No. represents the number of customers; D1-D3 represent demand of products 1-3.Figure 5
Logistics center and customer location in Chengdu.


### 6.1. Location Decision Results

The locations of the logistics center and customers (blue icon) are shown in Figure5. After the location decision, three regional logistics centers, eight city logistics centers, nine dealers may not operate. The non-operating logistics centers are identified by red rectangle. Reasons for not operating are analyzed.Two national logistics centers (numberedAi, yellow icon): Since the national distribution center needs to distribute to regional logistics centers and large customers, one national distribution center’s capacity is not enough, both of them need to operate.Four regional logistics centers (numberedBi, purple icon): Regional logistics center B2 needs to operate. B2 is closest to the two national logistics centers and is closer to the city logistics center compared with B1, B3, and B4. In the vicinity of B1, B3, and B4, the city logistics centers (numbered Ci) are less distributed and dispersed. Considering the operating costs, B1, B3, and B4 may not operate.Twelve city logistics centers (numberedCi, green icon): City logistics centers C1, C4, C8, and C11 need to operate. C2, C9, and C10 are far away from B2 and have high transportation costs compared with C4 and C8, so they are not operated. C3, C5, C6, C7, and C12 are remote and scattered, and the demand of dealer is small. It does not operate due to transportation and operating costs.Twenty-two dealers (numberedDi, red icon): Through the calculation of location and routing, 9 dealers do not need to operate. The remaining 13 dealers can fulfill the product distribution demands of 78 second-class customers.The above analyses show that M company's “2-4-12-22” supply chain network is too large. The current demand can be satisfied by the “2-1-4-13” scale network. The location-routing decision can reduce the operations of three regional distribution centers, eight city distribution centers, and nine dealers.


### 6.2. Transportation Route Decision Results

Distribution route of the three products for 85 customers is shown in Table4.Table 4
Delivery route of 3 products for 85 customers (Ton).
No. Delivery path No Delivery path No. Delivery path No. Delivery path 1 A1-B2-C1-D8-E1 23 A1-B2-C1-D7-E23 43 A1-B2-C11-D19-E43 67 A1-B2-C1-D1-E67 2 A1-B2-C1-D1-E2 24 A1-B2-C11-D2-E24 44 A1-B2-C1-D1-E44 68 A1-B2-C8-D21-E68 3 A1-B2-C1-D1-E3 25 A1-B2-C11-D2-E25 45 A1-B2-C11-D2-E45 69 A1-B2-C1-D1-E69 4 A1-B2-C11-D19-E4 26 A1-B2-C1-D1-E26 46 A1-B2-C1-D8-E46 70 A1-B2-C11-D19-E70 5 A1-B2-C1-D1-E5 27 A1-B2-C11-D2-E27 47 A1-B2-C11-D19-E47 71 A1-B2-C1-D1-E71 6 A1-B2-C11-D19-E6 28 A1-B2-C11-D19-E28 48 A1-B2-C11-D2-E48 72 A1-B2-C1-D1-E72 7 A1-B2-C1-D10-E7 29 A1-B2-C1-D8-E29 49 A1-B2-C1-D1-E49 73 A1-B2-C1-D1-E73 9 A1-B2-C1-D1-E9 30 A1-B2-C1-D1-E30 51 A1-B2-C8-D21-E51 74 A1-B2-C1-D1-E74 10 A1-B2-C1-D1-E10 31 A1-B2-C1-D1-E31 52 A1-B2-C1-D7-E52 75 A1-B2-C1-D1-E75 11 A1-B2-C1-D22-E11 32 A1-B2-C1-D1-E32 53 A1-B2-C1-D1-E53 76 A1-B2-C1-D1-E76 12 A1-B2-C1-D1-E12 33 A1-B2-C11-D2-E33 54 A1-B2-C1-D1-E54 77 A1-B2-C1-D1-E77 13 A1-B2-C8-D21-E13 34 A1-B2-C11-D19-E34 55 A1-B2-C1-D1-E55 78 A1-B2-C1-D1-E78 14 A1-B2-C1-D1-E14 35 A1-B2-C1-D22-E35 56 A1-B2-C1-D1-E56 79 A1-E79 15 A1-B2-C1-D1-E15 36 A1-B2-C1-D16-E36 57 A1-B2-C11-D17-E57 80 A2-E80 16 A1-B2-C11-D2-E16 37 A1-B2-C11-D2-E37 58 A1-B2-C1-D1-E58 81 A2-E81 17 A1-B2-C1-D22-E17 38 A1-B2-C11-D2-E38 60 A1-B2-C1-D8-E60 82 A2-E82 19 A1-B2-C1-D10-E19 39 A1-B2-C11-D19-E39 61 A1-B2-C1-D1-E61 83 A1-E83 20 A1-B2-C8-D21-E20 40 A1-B2-C4-D18-E40 62 A1-B2-C8-D21-E62 84 A2-E84 21 A1-B2-C8-D21-E21 41 A1-B2-C4-D18-E41 63 A1-B2-C1-D1-E63 85 A2-E85 22 A1-B2-C1-D8-E22 42 A1-B2-C1-D1-E42 64 A1-B2-C1-D4-E64 8 A1-B2-C1-D8-E8 8 A1-B2-C1-D8-E8 8 / 18 / 18 A1-B2-C1-D8-E18 18 A1-B2-C1-D8-E18 50 A1-B2-C8-D21-E50 50 A1-B2-C8-D21-E50 50 / 59 A1-B2-C4-D12-E59 59 A1-B2-C4-D12-E59 59 / 65 A1-B2-C1-D4-E65 65 A1-B2-C1-D4-E65 65 / 66 A1-B2-C4-D18-E66 66 A1-B2-C4-D18-E66 66 /If distribution routes of the three products are the same, the route decision follows the customer number, and the first 20 lines show the delivery of three products for 79 customers. The delivery routes of three products for 6 customers in the last 6 lines are not completely consistent and are listed from left to right in the order of delivery paths of products 1-3. This example shows that the mathematical model and the solution algorithm can complete both operational decision and routing decision of three products through five supply chain levels.


### 6.3. Sensitivity Analysis

In the case of 1%, 5%, and 10% change in demand, the impact of uncertain demand on total costs was analyzed. The total operating costs of the four demand changes under each experimental scale were counted. Demand change is represented byσ. Each row of Table 5 shows the total costs of different situations.Table 5
Comparison of costs and product volume when demand changes.
Scale of the study 1-1-2-4-10 2-2-4-10-35 2-2-6-14-55 2-4-12-22-85 3-5-18-70-205 (2+8) (4+31) (5+50) (7+78) (10+195) σ = 0 155118311.9 4861980932 5213955899 5363403216 5672995297 σ = 0.01 155478693.3 4913027119 5224634879 5350314683 5666793047 σ = 0.05 155446843.3 4913110137 5231902083 5392591433 5685517791 σ = 0.1 155867650.3 4941648089 5233955899 5389350467 5691250911 Mean 155477874.7 4907441569 5226112190 5373914950 5679139262 max-min 749338.4 79667157 20000000 42276750 24457864 (max-min)/mean 0.48% 1.62% 0.38% 0.79% 0.43%  Average value     0.74%As the demand changes, the total costs of the four scenarios change. As a result, the ratio of the maximum and minimum differences to the average cost is 0.74%. Therefore, considering each customer's demands for each product separately can reduce the impact of demand uncertainty on the total costs of supply chain network and reasonably predict the demand distribution.Based on the location decision results in Section6.1, another sensitivity analysis on location decision is conducted. Assuming that the model does not contain location decision, all dealers need to operate.The total costs and capacity without location decision are shown in Table6.Table 6
Comparison of non-site location and site location.
Total costs Total costs of customere1 Total costs of customere2 Product1 Product2 Product3 Product1 Product2 Product3 Time Volume of Volume of Volume of Volume of Volume of Volume of Customere1 Customere1 Customere1 Customere2 Customere2 Customere2 Non-site location 5367103216 3422230239 1944872977 2644 194 12 1815 427 192 68 Site location 5363403216 3422230239 1941172977 2644 194 12 1815 427 192 75 Difference 3700000 0 3700000 0 0 0 0 0 0 7The comparison of Table6 shows that the location decision can effectively save the total costs, which is 3,700,000 RMB, compared with the non-site selection decision. 3,700,000 RMB is exactly the sum of the operation costs of three regional logistics centers, eight city logistics, and nine dealers. Specific operational decisions are in Section 6.1.


## 6.1. Location Decision Results

The locations of the logistics center and customers (blue icon) are shown in Figure5. After the location decision, three regional logistics centers, eight city logistics centers, nine dealers may not operate. The non-operating logistics centers are identified by red rectangle. Reasons for not operating are analyzed.Two national logistics centers (numberedAi, yellow icon): Since the national distribution center needs to distribute to regional logistics centers and large customers, one national distribution center’s capacity is not enough, both of them need to operate.Four regional logistics centers (numberedBi, purple icon): Regional logistics center B2 needs to operate. B2 is closest to the two national logistics centers and is closer to the city logistics center compared with B1, B3, and B4. In the vicinity of B1, B3, and B4, the city logistics centers (numbered Ci) are less distributed and dispersed. Considering the operating costs, B1, B3, and B4 may not operate.Twelve city logistics centers (numberedCi, green icon): City logistics centers C1, C4, C8, and C11 need to operate. C2, C9, and C10 are far away from B2 and have high transportation costs compared with C4 and C8, so they are not operated. C3, C5, C6, C7, and C12 are remote and scattered, and the demand of dealer is small. It does not operate due to transportation and operating costs.Twenty-two dealers (numberedDi, red icon): Through the calculation of location and routing, 9 dealers do not need to operate. The remaining 13 dealers can fulfill the product distribution demands of 78 second-class customers.The above analyses show that M company's “2-4-12-22” supply chain network is too large. The current demand can be satisfied by the “2-1-4-13” scale network. The location-routing decision can reduce the operations of three regional distribution centers, eight city distribution centers, and nine dealers.


## 6.2. Transportation Route Decision Results

Distribution route of the three products for 85 customers is shown in Table4.Table 4
Delivery route of 3 products for 85 customers (Ton).
No. Delivery path No Delivery path No. Delivery path No. Delivery path 1 A1-B2-C1-D8-E1 23 A1-B2-C1-D7-E23 43 A1-B2-C11-D19-E43 67 A1-B2-C1-D1-E67 2 A1-B2-C1-D1-E2 24 A1-B2-C11-D2-E24 44 A1-B2-C1-D1-E44 68 A1-B2-C8-D21-E68 3 A1-B2-C1-D1-E3 25 A1-B2-C11-D2-E25 45 A1-B2-C11-D2-E45 69 A1-B2-C1-D1-E69 4 A1-B2-C11-D19-E4 26 A1-B2-C1-D1-E26 46 A1-B2-C1-D8-E46 70 A1-B2-C11-D19-E70 5 A1-B2-C1-D1-E5 27 A1-B2-C11-D2-E27 47 A1-B2-C11-D19-E47 71 A1-B2-C1-D1-E71 6 A1-B2-C11-D19-E6 28 A1-B2-C11-D19-E28 48 A1-B2-C11-D2-E48 72 A1-B2-C1-D1-E72 7 A1-B2-C1-D10-E7 29 A1-B2-C1-D8-E29 49 A1-B2-C1-D1-E49 73 A1-B2-C1-D1-E73 9 A1-B2-C1-D1-E9 30 A1-B2-C1-D1-E30 51 A1-B2-C8-D21-E51 74 A1-B2-C1-D1-E74 10 A1-B2-C1-D1-E10 31 A1-B2-C1-D1-E31 52 A1-B2-C1-D7-E52 75 A1-B2-C1-D1-E75 11 A1-B2-C1-D22-E11 32 A1-B2-C1-D1-E32 53 A1-B2-C1-D1-E53 76 A1-B2-C1-D1-E76 12 A1-B2-C1-D1-E12 33 A1-B2-C11-D2-E33 54 A1-B2-C1-D1-E54 77 A1-B2-C1-D1-E77 13 A1-B2-C8-D21-E13 34 A1-B2-C11-D19-E34 55 A1-B2-C1-D1-E55 78 A1-B2-C1-D1-E78 14 A1-B2-C1-D1-E14 35 A1-B2-C1-D22-E35 56 A1-B2-C1-D1-E56 79 A1-E79 15 A1-B2-C1-D1-E15 36 A1-B2-C1-D16-E36 57 A1-B2-C11-D17-E57 80 A2-E80 16 A1-B2-C11-D2-E16 37 A1-B2-C11-D2-E37 58 A1-B2-C1-D1-E58 81 A2-E81 17 A1-B2-C1-D22-E17 38 A1-B2-C11-D2-E38 60 A1-B2-C1-D8-E60 82 A2-E82 19 A1-B2-C1-D10-E19 39 A1-B2-C11-D19-E39 61 A1-B2-C1-D1-E61 83 A1-E83 20 A1-B2-C8-D21-E20 40 A1-B2-C4-D18-E40 62 A1-B2-C8-D21-E62 84 A2-E84 21 A1-B2-C8-D21-E21 41 A1-B2-C4-D18-E41 63 A1-B2-C1-D1-E63 85 A2-E85 22 A1-B2-C1-D8-E22 42 A1-B2-C1-D1-E42 64 A1-B2-C1-D4-E64 8 A1-B2-C1-D8-E8 8 A1-B2-C1-D8-E8 8 / 18 / 18 A1-B2-C1-D8-E18 18 A1-B2-C1-D8-E18 50 A1-B2-C8-D21-E50 50 A1-B2-C8-D21-E50 50 / 59 A1-B2-C4-D12-E59 59 A1-B2-C4-D12-E59 59 / 65 A1-B2-C1-D4-E65 65 A1-B2-C1-D4-E65 65 / 66 A1-B2-C4-D18-E66 66 A1-B2-C4-D18-E66 66 /If distribution routes of the three products are the same, the route decision follows the customer number, and the first 20 lines show the delivery of three products for 79 customers. The delivery routes of three products for 6 customers in the last 6 lines are not completely consistent and are listed from left to right in the order of delivery paths of products 1-3. This example shows that the mathematical model and the solution algorithm can complete both operational decision and routing decision of three products through five supply chain levels.


## 6.3. Sensitivity Analysis

In the case of 1%, 5%, and 10% change in demand, the impact of uncertain demand on total costs was analyzed. The total operating costs of the four demand changes under each experimental scale were counted. Demand change is represented byσ. Each row of Table 5 shows the total costs of different situations.Table 5
Comparison of costs and product volume when demand changes.
Scale of the study 1-1-2-4-10 2-2-4-10-35 2-2-6-14-55 2-4-12-22-85 3-5-18-70-205 (2+8) (4+31) (5+50) (7+78) (10+195) σ = 0 155118311.9 4861980932 5213955899 5363403216 5672995297 σ = 0.01 155478693.3 4913027119 5224634879 5350314683 5666793047 σ = 0.05 155446843.3 4913110137 5231902083 5392591433 5685517791 σ = 0.1 155867650.3 4941648089 5233955899 5389350467 5691250911 Mean 155477874.7 4907441569 5226112190 5373914950 5679139262 max-min 749338.4 79667157 20000000 42276750 24457864 (max-min)/mean 0.48% 1.62% 0.38% 0.79% 0.43%  Average value     0.74%As the demand changes, the total costs of the four scenarios change. As a result, the ratio of the maximum and minimum differences to the average cost is 0.74%. Therefore, considering each customer's demands for each product separately can reduce the impact of demand uncertainty on the total costs of supply chain network and reasonably predict the demand distribution.Based on the location decision results in Section6.1, another sensitivity analysis on location decision is conducted. Assuming that the model does not contain location decision, all dealers need to operate.The total costs and capacity without location decision are shown in Table6.Table 6
Comparison of non-site location and site location.
Total costs Total costs of customere1 Total costs of customere2 Product1 Product2 Product3 Product1 Product2 Product3 Time Volume of Volume of Volume of Volume of Volume of Volume of Customere1 Customere1 Customere1 Customere2 Customere2 Customere2 Non-site location 5367103216 3422230239 1944872977 2644 194 12 1815 427 192 68 Site location 5363403216 3422230239 1941172977 2644 194 12 1815 427 192 75 Difference 3700000 0 3700000 0 0 0 0 0 0 7The comparison of Table6 shows that the location decision can effectively save the total costs, which is 3,700,000 RMB, compared with the non-site selection decision. 3,700,000 RMB is exactly the sum of the operation costs of three regional logistics centers, eight city logistics, and nine dealers. Specific operational decisions are in Section 6.1.


## 7. Conclusion

This paper studies the joint optimization of location and routing problem with the background of distribution under Omni-Channel integration. In our research, a many-to-many distribution of supply chain network is designed. We consider each customer's uncertain demands for different products, then build a demand network, and complete the integration optimization of distribution network and demand network. The progress and contributions of this research include the following:( 1 ) The customers are classified in this designed supply chain network. This study increases product categories and increases supply chain levels. Many-to-many distribution is achieved between the distribution network levels, and the demand network considers the different demands of customers for different products. Less attention is concentrated on integration optimization issues of distribution network and demand network in existing research.( 2 ) Through three improvements of particle swarm optimization, such as collaborative dimension reduction, boundary processing, and introduction of mutation factors, the performance of the algorithm is improved. The paper provides an algorithmic solution to solve high-dimensional problems.In future study, the research results of this paper can be extended to the scale economy of the distribution network, while increasing the uncertainty factors and solving larger scale examples.



---

*Source: 1520213-2019-04-04.xml*