Cooperative Coevolutionary Bare-Bones Particle Swarm Optimization With Function Independent Decomposition for Large-Scale Supply Chain Network Design With Uncertainties

Supply chain network design (SCND) is a complicated constrained optimization problem that plays a significant role in the business management. This article extends the SCND model to a large-scale SCND with uncertainties (LUSCND), which is more practical but also more challenging. However, it is difficult for traditional approaches to obtain the feasible solutions in the large-scale search space within the limited time. This article proposes a cooperative coevolutionary bare-bones particle swarm optimization (CCBBPSO) with function independent decomposition (FID), called CCBBPSO-FID, for a multiperiod three-echelon LUSCND problem. For the large-scale issue, binary encoding of the original model is converted to integer encoding for dimensionality reduction, and a novel FID is designed to efficiently decompose the problem. For obtaining the feasible solutions, two repair methods are designed to repair the infeasible solutions that appear frequently in the LUSCND problem. A step translation method is proposed to deal with the variables out of bounds, and a labeled reposition operator with adaptive probabilities is designed to repair the infeasible solutions that violate the constraints. Experiments are conducted on 405 instances with three different scales. The results show that CCBBPSO-FID has an evident superiority over contestant algorithms.


I. INTRODUCTION
T HE SUPPLY CHAIN network (SCN) has attracted increasing attention in different fields, such as production manufacturing [1], [2] and business management [3], [4].The members of SCN include suppliers, manufacturers, warehouses, distributors, retailers, and customers.These members and the flows of finances and materials among them make up the entire network.The aim of the design of SCN (SCND) is to minimize the total cost of the entire system and to meet demands of different members [5]- [7].An excellent SCND can optimize enterprise management, improve efficiency, and finally encourage consumptions.
Different SCND models have different features.For example, SCND under uncertainties (USCND) has uncertain factors, such as uncertain demand of customers and uncertain supply of suppliers [8].These uncertain factors make it difficult to evaluate the solutions and generate feasible solutions.There are two main simulation methods of uncertainties to help the fitness evaluation (FE): 1) the scenario analysis and 2) the Monte Carlo method (MCM).The scenario analysis only analyzes limited scenarios of uncertainties [9].Differently, MCM simulates the uncertainties in each evaluation of solutions and is more complicated [10].This article also adopts MCM to the USCND.
The SCND/USCND problem can be regarded as a constrained optimization problem (COP).Due to the great potential of evolutionary algorithms (EAs) in solving COPs, some researchers have already used EAs for SCND problems [11].Salem and Haouari [10] applied the particle swarm optimization (PSO) to solve a more complex USCND problem with multiple periods and uncertainties in both supply and demand.Sahebjamnia et al. [12] used the red deer algorithm (RDA) for the sustainable tire closed-loop SCND problem.
However, it is still difficult for the existing algorithms to effectively solve the large-scale USCND (LUSCND) problem, which involves both the large-scale challenge and the uncertain challenge.On the large-scale issue, a large number of suppliers and customers are included in the LUSCND to meet the growing market demand [13].Unfortunately, these increasing SCN members lead to a larger solution space and more constraints, making it more difficult to find feasible and optimal solutions.An efficient method for the large-scale issue is in great need to obtain the satisfactory solutions within the limited time.Besides, as the scale of the LUSCND becomes large, the number of constraints will increase sharply, and infeasible solutions may appear frequently in this complex problem.Most of the current works for SCND only apply penalty functions on infeasible solutions, being useless to find a feasible solution for the LUSCND problem.Moreover, in the uncertain environment, the constraints in the LUSCND are also uncertain [10].These uncertain constraints will also increase the number of constraints.Therefore, it is hard for traditional approaches to obtain feasible solutions of the LUSCND problem.New methods are needed to replace the penalty functions and to help finding the feasible solutions.
Therefore, this article focuses on the practical and challenging LUSCND problem and proposes an efficient approach to solve this problem.On the aspect of the problem model, the USCND in [10] is extended with large-scale variables.On the aspect of the algorithm design, an efficient method for the large-scale issue and two new methods for constraint handling are proposed for the LUSCND problem with complex constraints.
First, to relieve the large-scale challenges, the model is simplified by replacing the binary coding of the location with the integer coding (under the assumption that a member, like a supplier, can choose one location at most), which helps to reduce the solution space of the problem.If there are n binary values (n locations to be chosen), the solution space of the location will reduce from 2 n to n.Moreover, a cooperative coevolutionary bare-bones PSO (CCBBPSO) with function independent decomposition (FID) is proposed to solve the LUSCND problem.The cooperative coevolutionary (CC) strategy is a promising solver for large-scale optimization problems [14]- [16].It decomposes the problems into several parts and solves these parts independently via different populations, which can search solutions in different directions and help to increase the solution diversity.Although there have been many decomposition methods, determining how to decompose problems effectively is still a critical problem for the CC strategy.Differential grouping (DG) via interdependencies between variables (IaV) has a high grouping accuracy [15], [16].However, the IaV is always difficult to be obtained because it needs complex calculations and consumes many FEs.On the contrary, the problem structure of LSCND is relatively specific because the variables have different functions.That is, the variables of LSCND belong to different kinds of SCN members which perform different tasks in the problem.Therefore, this article proposes an FID-based CC strategy to decompose the LSCND problem according to different functions of variables.
Second, to deal with the infeasible solutions of the LUSCND problem, two repair methods are designed.In the traditional works for COPs, the gradient-based repair method is effective to repair the infeasible solutions by the gradient information from the infeasible solutions to the feasible region [17], [18].However, it is extremely difficult to obtain the gradient information of the LUSCND problem because there are a large number of variables.In this article, a step translation method (STM) is designed to map the out-of-bound values within the boundaries, and a labeled reposition operator with adaptive probabilities (LRO ap ) is designed to repair infeasible solutions that violate the constraints.
The main contributions of this article are shown as follows.
1) A new encoding scheme for the LUSCND model is proposed to reduce the solution space.
2) The CCBBPSO with FID is proposed to solve the LUSCND problem.During the evolution, the problem is decomposed into several subproblems by FID, and these subproblems are solved independently.3) STM and LRO ap are designed to repair infeasible solutions and improve the solving efficiency of the CCBBPSO-FID algorithm.These innovative strategies enable CCBBPSO-FID to deal with larger data scales from thousands to hundreds of thousands dimensions.To evaluate the performance of CCBBPSO-FID, three scales of the instance sets are generated.The proposed algorithm is compared with several existing algorithms.First, as the LUSCND model is extended from the USCND model in [10], the global version PSO (GPSO) used in [10] is compared.Also, BBPSO [19] is compared because it is the basic algorithm of CCBBPSO-FID.Second, the competitive swarm optimizer (CSO) [20] is compared because it is an efficient algorithm for large-scale optimization.Moreover, two novel algorithms, including the social engineering optimizer (SEO) [21] and RDA [12], are compared because the SEO is a recent algorithm for complex problems and the RDA is a recent algorithm used for the LUSCND problems.Third, CCBBPSO-FID is compared with non-PSO (e.g., differential evolution (DE)-based [22]) algorithms and CCBBPSO with other decomposition schemes to verify the advantages by combining CCBBPSO with FID.For the fair comparison, these algorithms are all combined with the proposed repair methods, including STM and LRO ap .Fourth, all algorithms without the repair methods are compared to prove the effectiveness of STM and LRO ap .The experimental results show that CCBBPSO-FID obtains a better performance on the LUSCND problem.
The remainder of this article is organized as follows.Section II introduces the definitions of LUSCND and the original BBPSO algorithm.Section III describes all components of CCBBPSO-FID in detail.Section IV presents the exhaustive experiments and analyzes the results in different aspects.Finally, Section V gives a conclusion of this article.

A. LUSCND Problem Definition
SCND is to design a network which is composed of a logistic system, a supply system, and a distribution system [23].
The logistic system records the order information of customers and warehouses; the supply system reflects the production patterns of suppliers; and the distribution system describes how products are distributed among different members.This article focuses on the LUSCND with a large-scale logistic system under supply and demand uncertainties [10].There have been some relevant studies on the LUSCND problem, and the main similarities and differences between them and this article are shown in Table S.I in the supplementary material, from both the model formulation aspect and the algorithm design aspect.The main similarity is that the algorithms used in these studies are all EAs, and the main differences include different problem models and algorithm design details.As shown in Table S.I in the supplementary material, this article mainly focuses on and contributes to the algorithm design like the large-scale tackling and the constraint handling to propose an efficient algorithm but not the new problem model.In this sense, we adopt the USCND model in [10] because it considers multiple periods.Moreover, this model considers only the economic dimension so that it fits our objective to minimize the total cost of the entire system.Nevertheless, we extend the USCND to the large-scale model to make it more practical.To make the LUSCND model more reasonable, the following suppositions and constraints are considered in this article.
1) Different from most of the existing works that consider only one period, our model considers multiple periods and three echelons, including suppliers, warehouses, and customers.2) This model has the capacity limitation.For example, the order quantity of a supplier cannot exceed its capacity; the inventory and the order quantity of a warehouse cannot exceed its capacity; and the supply of a customer cannot exceed its demand.3) A supplier or warehouse can choose one location at most.4) Each supplier can supply for several warehouses, and each warehouse can transport to several customers, and vice-versa.This model is applicable in the three-echelon LUSCND that includes suppliers, warehouses, and customers.Fig. 1 illustrates the application scenario of this LUSCND with S suppliers, W warehouses, and C customers.In this model, the workflow mainly includes four stages.
1) Customers send the uncertain demands to warehouses (C_demand) in a continuous period of time.2) Some warehouses with the optimal locations (W_open) are chosen to serve the customers and will send order to the suppliers (WS_order).3) Some suppliers with the optimal locations (S_open) are chosen to serve warehouses and provide products to warehouses (WS_order × S_yield, where S_yield represents the uncertain output rates of suppliers).4) Warehouses chosen in the second stage will transport the products to the customers (WC_ship).It should be noted that as shown in Fig. 1, the uncertainties include the demand of customers and the supply of suppliers to make the model closer to the practical applications.To satisfy the uncertain demand of customers, it is better to supply materials as much as possible, but too much supply will cost more, including the fixed cost, the inventory cost, and the shipping cost.Therefore, to satisfy the demand and cut the cost simultaneously, the competent locations of suppliers and warehouses (S_open and W_open) should be optimally chosen.Moreover, the order quantity of warehouses to suppliers (WS_order) and the transport quantity of warehouses to customers (WC_ship) should be reasonably optimized to satisfy the demand of customers with the minimum cost.
In order to clearly describe the relationship between the decision variables and the objective function of the model, the details of the model formulation are described as follows.For customer k(k = 1, 2, . . ., C), in period t(t = 1, 2, . . ., T), it has a random demand quantity C_demand k,t and the penalty C_penalty k for a unit of the demand unmet.For warehouse j(j = 1, 2, . . ., W), it has W_n locations that can be chosen, and the lth (l = 1, 2, . . ., W_n) location is with the capacity of W_capacity j,l and the fixed cost W_fixedCost j,l (e.g., the rent cost, the maintenance cost, etc.).The inventory cost is W_inventoryCost j for storing a unit of goods in a period, and the initial inventory is W_inventory j,0 .For supplier i (i = 1, 2, . . ., S), it has S_n locations to be chosen.The capacity of the rth (r = 1, 2, . . ., S_n) location is S_capacity i,r and the fixed cost is S_fixedCost i,r (e.g., the production cost, the management cost, etc.).Besides, if a unit of capacity of supplier i is unused in a period, the penalty is S_penalty i .The random output rate of supplier i at period t is denoted by S_yield i,t , and the real supply quantity is the total order quantity of warehouses multiplied by S_yield i,t .The transportation costs from supplier i to warehouse j and from warehouse j to customer k with a unit of goods are denoted by T_sup2war i,j and T_war2cus j,k , respectively.
The decision variables include S_open (the chosen locations of suppliers), W_open (the chosen locations of warehouses), WS_order (the order quantity of warehouses to suppliers), and WC_ship (the shipping quantity of warehouses to customers).S_open i,r represents whether the rth location of supplier i is selected (S_open i,r = 1) or not (S_open i,r = 0), and W_open j,l represents whether the lth location of warehouse j is selected (W_open j,l = 1) or not (W_open j,l = 0).WS_order i,j,t represents the order quantity of supplier i from warehouse j at period t, and WC_ship j,k,t represents the shipping quantity of warehouse j to customer k at period t.The objective is to obtain the optimal S_open, W_open, WS_order, and WC_ship to minimize the total cost.The structure of the solution is shown in Fig. 2, where S_open and W_open are binary numbers, and WS_order and WC_ship are real numbers.
The total cost totalCost of LUSCND can be divided into four parts, including the fixed cost, the inventory cost, the shipping cost, and the penalty.Specifically, the fixed cost fixedCost can be formulated as The inventory cost inventoryCost can be formulated as where W_inventory j,t is the inventory quantity of warehouse j at period t, and it is calculated as follows to keep the production-inventory balance of warehouse j: where S_output i,j,t is the real output quantity of supplier i to warehouse j at period t and is calculated as follows: The shipping cost shippingCost can be formulated as T_war2cus j,k × WC_ship j,k,t . ( At last, the penalty Penalty can be formulated as where S_rest i,t is the unused capacity of supplier i at period t, and C_unmet k,t is the unmet quantity of customer k at period t, and they can be calculated as follows: Therefore, the final objective function of LUSCND can be formulated as Minimize totalCost = fixedCost + inventoryCost subject to: The model formulation is clarified according to the workflow of the LUSCND problem, as shown in Fig. 1.In the first stage, customers send demands to warehouses (C_demand) in a continuous period of time (T).In the second stage, some warehouses with the optimal locations are chosen to serve the customers (W_open) and will send order to suppliers (WS_order).In constraint (11), if W_n l=1 W_open j,l = 0, warehouse j is not open; otherwise, the lth location with W_open j,l = 1 represents the location that warehouse j selects.The order quantity of warehouses to suppliers (WS_order) should not be smaller than 0, as shown in (12), and should be not larger than the capacity of suppliers (S_capacity), as shown in (14).In addition, the current inventory of warehouses [W_inventory calculated by ( 3)] and the order quantity of warehouses to suppliers (WS_order) should not be larger than the capacity of warehouses (W_capacity), as shown in (15).The current inventory of warehouses (W_inventory) should not be smaller than 0, as shown in (16).In the third stage, some suppliers with the optimal locations are chosen to serve warehouses (S_open) and provide products to warehouses (S_yield × WS_order).In constraint (10), if S_n r=1 S_open i,r = 0, supplier i is not open; otherwise, the rth location with S_open i,r = 1 represents the location that supplier i selects.In the fourth stage, the warehouses chosen in the second stage will transport the products to the customers (WC_ship).The transport quantity (WC_ship) should not be smaller than 0, as shown in (13).In addition, the unmet quantity of customers (C_demand -WC_ship) should not be smaller than 0, as shown in (17).
Table I summarizes the notations of all the indices, parameters, and decision variables used in the LUSCND.To make the real-world applications of our proposed approach clearer for The example can be stated in four steps.First, customers send demands to warehouses.It is assumed that the total demand quantities of the four customers are 3, 2, 4, and 3, respectively.It is noted that there is no need to calculate the demand quantity of each customer to each warehouse in the FE, since the related cost (the penalty from customers) can be obtained by WC_ship in ( 6) and (8).Second, warehouses with the optimal locations are chosen to serve the customers (W_open, W 1,3 and W 2,2 in this example), and send orders to the supplier (WS_order).For example, W 1,3 sends 2 units of order to S 1,3 .Third, the supplier with the optimal location is chosen (S_open, S 1,3 in the example) and transports goods to warehouses (S_yield × WS_order), shown as the dotted lines with a number on the right.The yield rate of S 1,3 is 0.8, and it can deliver 2 × 0.8 = 1.6 units of goods to W 1,3 .Similarly, S 1,3 sends 3 × 0.8 = 2.4 units of goods to W 2,2 .After getting the supply, the quantities of W 1,3 and W 2,2 including supplies and inventory are 1.6 + 3 = 4.6 and 2.4 + 2 = 4.4, respectively.Finally, the warehouses that have been chosen in the second step deliver goods to customers (WC_ship).As shown in the parentheses, C 1 , C 2 , C 3 , and C 4 get 3, 2, 2.6, 1.4 units of goods, respectively.The demands of C 3 and C 4 are not satisfied.

B. BBPSO
BBPSO, developed by Kennedy [19], is a simplified version of PSO.Instead of using the velocity vector to update the position vector, BBPSO applies the Gaussian distribution to set the position information directly, as shown in (18) x i,d (g where the Gaussian distribution with the mean value μ i,d (g) and the deviation σ i,d (g); pbest i,d (g) is the personal best of the ith individual; and gbest d (g) is the dth dimension of the global best individual.Another advantage of BBPSO is that there are no parameters needed to be tuned.Thus, BBPSO has attracted great attention in diverse applications [24]- [27].

A. Solution Encoding for Dimensionality Reduction
According to the LUSCND problem definition, the dimension of a solution (including S_open, W_open, WS_order, and WC_ship) is equal to When the LUSCND problem consists of tens of suppliers, warehouses, and periods and one hundred customers, the dimension increases sharply to 10 000.It needs a large storage for a solution, which increases the problem-solving difficulty.An effective approach is to reduce the dimensionality of the model.From the view of data types of decision variables, WS_order and WC_ship with real numbers cannot be compressed.For S_open ( S_n r=1 S_open i,r ∈ {0, 1}) with binary numbers, the array [S_open i,1 , S_open i,2 , . . ., S_open i,S_n ] ∈ {0, 1} S_n can be replaced by an integer S_open i ∈ {0, 1, 2, . . ., S_n}.If S_open i is 0, supplier i is not open; otherwise, S_open i represents the location that supplier i selects.Similar to S_open, [W_open j,1 , W_open j,2 , . . ., W_open j,W_n ] ∈ {0, 1} W_n can be replaced by W_open j ∈ {0, 1, 2, . . ., W_n}.Then, the dimension will be reduced from S × S_n + W × W_n to S + W. Meanwhile, the runtime will be reduced during the entire evolutionary process.Fig. 4 illustrates an example of the encoding scheme of S_open and WS_order in this article.Because the real number coding is used in the programming, values in S_open are needed to be rounded down to be feasible.For example, in Fig. 4, after being rounded down, S_open 1 = 3 represents that the first supplier (S 1 ) is open in the third location, and S_open 2 = 0 represents that the second supplier (S 2 ) is not open.Values in WS_order represent the order quantity of warehouses to suppliers.For example, WS_order 1,2,0 = 3.6 (in the first row and the second column) represents that the first warehouse (W 1 ) sends an order of 3.6 units to the second supplier (S 2 ) at the period t = 0.The encoding schemes of W_open and WC_ship are similar to these of S_open and WS_order, respectively, and they are not shown in the figure.After these modifications and removing the intermediate variables, the final model formulation of the objective function is more clearly shown in (21), as shown at the top of the previous page, whose relationship with the decision variables, including S_open, W_open, WS_order, and WC_ship, is also clearer.Moreover, the objective function is subject to ( 12), (13), and ( 22)- (27), as shown at the top of the previous page.Constraints ( 22)-( 27) are equivalent to (10), (11), and ( 14)- (17), respectively.In order to minimize the totalCost of ( 21), the locations of suppliers (S_open) have to be optimally chosen so that they not only satisfy ( 22) and ( 24), but also have low S_fixedCost and small S_capacity.Also, the locations of warehouses (W_open) have to be optimally chosen so that they satisfy ( 23) and (25) and have low W_fixedCost.Moreover, the WS_order and WC_ship have to be optimized to be subject to (12), (13), and ( 24)- (27).

B. Fitness Evaluation With Uncertainties
Apart from the large-scale challenge in the LUSCND, how to deal with the uncertainties is an another difficulty.The FEs of solutions in the uncertain problems and the dynamic optimization problems (DOPs) are both not fixed during the entire evolution of algorithms.It is noted that the fitness function of DOPs in each certain environment is fixed although it changes in different environments [28]- [30].However, in the LUSCND problem, the fitness function is always uncertain, with the uncertain demand of customers (C_demand) and the uncertain supply of suppliers (S_yield).Therefore, the strategies for DOPs are not suitable for the LUSCND problem.In this article, the MCM is used to evaluate the LUSCND problem by multiple samplings of the uncertain factors.
In the literature, the MCM, also called statistic testing method or random sampling, is used to solve mathematical and physical problems with uncertainties [31]- [33].It simulates uncertain factors through random sampling and obtains approximate fitness values.Due to its simplicity, this method is widely applied in the USCND problems [10], [34].Therefore, the MCM is also used to deal with the supply and demand uncertainties in this article.For each solution, ten simulations with different S_yield and C_demand are randomly sampled.The fitness value of each simulation is calculated as (21) by using certain values of uncertain factors.Then, the mean result of the ten simulations is regarded as the final fitness value.
It is worth mentioning that fixedCost and the shipping cost of warehouses to customers (shippingCost2) are not related to uncertain factors, so there is no need to recalculate them in the simulation.Other costs, including inventoryCost, the shipping cost of suppliers to warehouses (shippingCost1), and Penalty, are calculated ten times with ten different samples of S_yield and C_demand.The sum of fixedCost and shipping-Cost2, and the mean result of ten simulations is recorded as the final fitness value.During the calculation, the feasibility of the solution is judged, and an error code (ecode) is recorded which will be used in the repair operator LRO ap in Section III-F.If ecode is equal to 0, it means that the individual is feasible, and the calculation continues until the MCM is executed ten times; otherwise, it means that the individual is infeasible, and this calculation will end in advance.The fitness value of infeasible solutions is set to 10 000 000. Fig. 5 shows two cases in the calculation, including a feasible solution and an infeasible solution.It is noted that Fig. 5(b) illustrates an infeasible solution case, where the number 7 is only an example and it can also be any other values smaller than 10.

C. Initialization
Random initialization is used to increase the population diversity and distribute individuals among the solution space widely.For S_open and W_open, they are generated randomly within {0, 1, 2, . . ., S_n} and {0, 1, 2, . . ., W_n}, respectively.For WS_order and WC_ship, two generation methods are designed, and the common point of them is to ensure that WS_order and WC_ship are smaller than S_capacity and W_capacity, respectively.The first method is that they are set randomly within [0, capacity], where capacity is a known parameter to set S_capacity and W_capacity.The second method is to deduce their values from formulas, being more complicated.From (24), it is deduced that WS_order can be set randomly within [0, capacity/W].In (3) and ( 4), if S_yield i,t is 1 and the inventory of warehouse j (including W_inventory j,t and W_inventory j,t−1 ) is 0, S i=1 WS_order i,j,t ≈ C k=1 WC_ship j,k,t and WC_ship can be set randomly within [0, (capacity/W × S)/C].
The initialization includes two steps.First, the global optimal solution xbest is generated randomly by the first method.Second, four subpopulations (swarm 1 , swarm 2 , swarm 3 , and swarm 4 ) are set separately via the second method according to the construction method introduced in Section III-D.For example, a solution x in swarm 1 consists of x.S_open (initialized in swarm 1 ) and xbest.W_open, xbest.WS_order, xbest.WC_ship.

D. Cooperative Coevolution
Inspired by the idea of "divide-and-conquer," Potter and De Jong [35] designed the CC strategy.It is an effective method for large-scale optimization [14]- [16], [36]- [38].The first step of this strategy is to divide the variables into several parts with a specific decomposition strategy, which can be regarded as the dimensionality reduction.Then, these parts are solved with different subpopulations independently, and the local best solutions of subpopulations are obtained to update xbest.
1) FID-Based Decomposition: Since the CC strategy was proposed, there have been several decomposition methods [14]- [16].The aim is to assign interacting variables to one group as more as possible.If interdependent variables are placed in different groups, they will be blind to the synchronized information of other interdependent variables, since different groups are evolved separately, and it will also be difficult to find the best values of these interdependent variables.The popular approach is to group variables via IaV [14]- [16], but it has high computational complexity, especially for super-large-scale problems.However, with the clear solution structure of LUSCND, there is no need to calculate IaV.
In this article, FID is devised to divide the problem according to different functions of different parts of the solution structure.The solution structure consists of four parts (S_open, W_open, WS_order, and WC_ship).Different parts have different functions, since they belong to different kinds of SCN members which perform different tasks.Specifically, the decision variables of S_open belong to suppliers which decide the locations of suppliers; the decision variables of W_open belong to warehouses which decide the locations of warehouses; the decision variables of WS_order belong to warehouses which decide the order quantities of warehouses to suppliers; and the decision variables of WC_ship belong to warehouses which decide the transportation quantities of warehouses to customers.Therefore, this article decomposes the LUSCND problem into four parts to deal with different subproblems, including locations of suppliers, locations of warehouses, order quantities of warehouses to suppliers, and transportation quantities of warehouses to customers.The four decomposed parts include S_open, W_open, WS_order, and WC_ship.In the proposed CCBBPSO algorithm, four different subpopulations, called swarm 1 , swarm 2 , swarm 3 , and swarm 4 , are used to evolve the four different parts, respectively.
Different from other decomposition methods that have to learn the IaV or to regroup variables during the evolution, FID based on different functions divides the problem naturally.It is executed only at the beginning of the evolution and does not consume FEs.In addition, FID has more general applicability in large-scale realistic problems, since it decomposes the solution by functions of different roles in problems.
2) Construction of Solutions in Subpopulations: Although decomposed parts of the solution space are evolved separately, it is necessary to constitute the entire solution to obtain the fitness value.In the proposed algorithm, a partial solution of an individual is integrated with other parts of xbest to obtain a complete solution, same as [14].For example, a solution x in swarm 1 is composed of x.S_open

E. STM-Based Boundary Processing
After updating the individual x by (18), if some variables exceed the search ranges, they will be restricted within reasonable ranges to satisfy the constraints.For variables belonging to S_open and W_open, they are limited within {0, 1, . . ., S_n} and {0, 1, . . ., W_n}, respectively.
For variables belonging to WS_order, such as WS_order i,j,t , it is limited within [0, min(S_capacity i,r , W_capacity j,l )], where i and j are the corresponding supplier and warehouse, and r and l are the corresponding locations.It indicates that the order quantity of supplier i from warehouse j at period t cannot exceed the capacity of warehouse j and supplier i.
For variables belonging to WC_ship, such as WC_ship j,k,t , it is limited within [0, min(W_capacity j,l , C_demand k,t )], where k and t are the corresponding customer and period.It shows that the shipping quantity of warehouse j to customer k at period t cannot exceed the capacity of warehouse j and the demand of customer k.
For variables out of bounds, if they are simply set to the bound values, many infeasible variables will have the same value (i.e., the bound value).This method also decreases the diversity of solutions.Therefore, STM is proposed in this article to map illegal values to legal values.The illustration of STM is shown in Fig. 6, where X is the illegal value and X is the corresponding legal value after the STM operation.Herein, L and U are the lower bound and the upper bound, respectively, and len is the length between L and U.Then, X is translated to X by the step of len, where − → XX is the motion vector.That is, if X > U, X will be moved to the left until it becomes smaller than U(X = X − ceil((X − U)/(U − L)) × (U − L)); otherwise if X < L, X will be moved to the right until it becomes larger than L(X = X + ceil((L − X)/(U − L)) × (U − L)).

F. LRO-Based Infeasible Solution Repairing
The penalty function is widely used in the complicated problems to calculate the fitness values of infeasible solutions [10], [39].However, with the scale growing, the number of constraints will rise sharply, making it more difficult to obtain a feasible solution.Therefore, the penalization is not an effective method, and the efficient way is to repair infeasible solutions.In this article, LRO is proposed to repair illegal values repeatedly by resetting them to reasonable values.LRO includes four repair operators Repair 1 , Repair 2 , Repair 3 , and Repair 4 , which are designed for constraints ( 24)-( 27), respectively.Only one kind of the operators is executed in each repair, and the label is set to record which repair operator will be executed.Recording the label helps to save running time, since the algorithm does not need to rejudge which repair operator to be used every time.For a complex problem with multiple constraints, LRO is an effective repair method to deal with these constraints separately, and the label helps to distinguish different repair operators for different constraints in the multiple repair process.
For simplicity, each repair operator only changes a specified part of solutions.Specifically, Repair 1 only repairs S_open; Repair 2 only repairs W_open; and Repair 3 and Repair 4 only repair WC_ship.Since WS_order is involved in ( 24)-( 26) and hard to update, it will not be changed in LRO.
In Repair 1 for ( 24), if W j=1 WS_order i,j,t > max r∈{1,...,S_n} (S_capacity i,r ), it represents the total order quantity of supplier i at period t exceeds the maximum capacity of supplier i, and this repair operator will be useless after changing S_open; otherwise, supplier i will choose one location where it can hold the capacity of W j=1 WS_order i,j,t .Considering the corresponding costs, the legal location with the smallest capacity and fixed cost is chosen.
In Repair 2 for (25), operations are similar to those in Repair 1 .If W_inventory j,t−1 + S i=1 WS_order i,j,t > max l∈{1,...,W_n} (W_capacity j,l ), where W_inventory j,t−1 is calculated by (3), it represents the sum of the previous inventory and new supplies of warehouse j at period t exceed the warehouse's maximum capacity, and this operator will not work after changing W_open; otherwise, warehouse j will choose the legal location with the smallest capacity and fixed cost.
If W_inventory j,t is smaller than 0, it represents constraint (16) [equivalent to (26)] is not satisfied, and Repair 3 will be used.Let difference denote −W_inventory j,t of the infeasible solution; then, W_inventory j,t should increase to 0 at least to make the solution feasible.Let WC_sum denote C k=1 WC_ship j,k,t .Since only WC_ship is changed in this operator, it can be deduced from (3) that WC_sum should decrease by difference.If WC_sum < difference, the operator is useless after updating WC_ship; otherwise, WC_ship j,k,t can decrease proportionally by difference × (WC_ship j,k,t /WC_sum) to satisfy the condition.It is remarkable that the precision problem also makes the theoretical feasible solutions infeasible.To solve the problem, WC_ship j,k,t should be subtracted a small real number precision (set as 0.02 in the algorithm), if WC_ship j,k,t is bigger than precision.Operations in Repair 4 for (17) [equivalent to (27)] are similar to those in Repair 3 .If C_unmet k,t is smaller than 0, it represents constraint ( 17) is not satisfied.Let difference denote −C_unmet k,t and WC_sum denote W j=1 WC_ship j,k,t in (8).Then, if WC_sum = difference, WC_ship j,k,t can decrease proportionally by difference × (WC_ship j,k,t /WC_sum).
After the FE, if ecode (the label) is e(e = 1, 2, 3, 4), Repair e will be executed.At the end of a repair operator, the fitness value is calculated again, and ecode is also returned.If ecode is 0, it represents that the solution turns to be feasible, and this repair operator ends with success; otherwise, if the solution is still infeasible after being repaired T (the number of periods in the LUSCND problem) times, LRO ends with failure.Fig. 7 shows two examples of LRO, including successful and unsuccessful repair.
LRO is useful but time consuming.In the proposed algorithm, P_repair i is set adaptively to control the repair probability of individual i, and the initial value is set to 0.01 × S (the number of suppliers in the LUSCND problem).If individual i has been repaired in the last generation, P_repair i will decrease by 0.001 × S; otherwise, this value will increase by 0.001 × S. LRO with the adaptive probability is called as LRO ap , and it helps to find feasible solutions of the largescale problem within the acceptable time, being proved in the experiments.

G. Complete CCBBPSO-FID Algorithm
The flowchart of CCBBPSO-FID is shown in Fig. 8.As mentioned above, the important parts of CCBBPSO-FID include the CC strategy based on FID, the boundary processing method STM, and the repair operator LRO ap .

A. Experimental Setting
It is difficult to obtain the instance set of the LUSCND problem in real life, especially for large-scale data, so the instance set is generated randomly.First, all capacities (S_capacity and W_capacity) are generated randomly within [capacity, 2 × capacity], where capacity is set to 5.0; all costs (S_penalty, S_fixedCost, W_inventoryCost, W_fixedCost, C_penalty, T_sup2war, and T_ware2cus) are set randomly within [0, cost], where cost is set to 10.0.Second, W_inventory j,0 (j = 1, 2, . . ., W), S_yield, and C_demand are set according to the following methods.
1) Generation of W_inventory j,0 : By analyzing ( 25) and ( 26), it can be observed that the feasibility of solutions is closely related to W_inventory j,0 .Equations ( 24) and ( 25) are considered to obtain reasonable values of W_inventory j,0 , as shown in Fig. 9.

B. Experimental Setup
All algorithms were implemented in C++, and experiments were run on the PC with Intel Core i7-7700 and 8.0-GB memory.Compared algorithms are classified into three groups: the first one includes PSOs, such as GPSO [10] and BBPSO [19]; the second one includes a largescale optimization algorithm-CSO [20] and two recent algorithms-SEO [21] and RDA [12].Furthermore, to verify the effectiveness of combining CCBBPSO with FID, the third group includes CCDE-FID (using DE/rand/1 in [22] as the optimizer) and two other CCBBPSO variants called average decomposition of CCBBPSO (CCBBPSO-AD) and random decomposition of CCBBPSO (CCBBPSO-RD).It is noted that four versions of DE in [22] are tested, including DE/rand/1 with the crossover rate CR = 0.1 and 0.9, and DE/best/1 with CR = 0.1 and 0.9 (the amplification factor F = 0.5 + 0.5 × rand(0, 1) for all versions).Herein, only the results of the first version (DE/rand/1 with CR = 0.1) are present because it has generally the best performance among the four DE variants.Besides, DG is not compared with FID since it consumes a large number of FEs, as shown in Table S.II in the supplementary material.STM and LRO ap are used in all algorithms for the fair comparison.
The control variable method [40] is used to investigate the parameters of all the algorithms to find out the optimal values,  and the results are shown as Tables S.III-S.V in the supplementary material.The optimal parameter settings of algorithms are shown in Table III.The settings of CCBBPSO-AD and CCBBPSO-RD are the same as that of CCBBPSO-FID.Each instance is tested 30 times independently, and the mean results are recorded.

C. Experimental Comparisons With Other Algorithms
Exhaustive experiments are conducted on 405 instances, including Test I, Test II, and Test III which are shown in Table II.Each instance set contains 135 (5 × 27) instances.For saving space, these instances (Inst.)are classified in three different ways: 1) based on the level of S_yield, each test set is divided into YL, YM, and YH; 2) based on the CV level of C_demand, each test set is divided into DL, DM, and DH; and 3) based on the PDF of C_demand, each test set is divided into Norm, Log, and Tri.Each group includes 45 (135/3) instances (each instance has an average result of 30 times), and the average result [totalCost in (21)] of each group is shown in Tables IV-VI.For example, to calculate the average result of the group "YL," we first calculate the average result of 30 times for each instance in YL, and then calculate the average result of 45 instances.The last row (Avg.) in tables shows the average result of each set.Results in italic type mean that the global best solution is infeasible in some tests.Results in bold type are the best results among all algorithms.Besides, all results of CCBBPSO-RD and the results of SEO and RDA on Test III are equal to 10 000 000 and, therefore, are not shown in tables, indicating that CCBBPSO-RD, SEO, and RDA cannot obtain feasible solutions on the corresponding test sets.IV, it can be seen that CCBBPSO-FID yields the minimum cost in each group with an obvious advantage, followed by CSO and SEO.On the contrary, CCBBPSO-AD and CCDE-FID obtain infeasible solutions in most of the test groups.As S in this dataset is small, the frequency of LRO ap (related to repair probability) drops, so some algorithms cannot obtain feasible solutions.GPSO, BBPSO, and RDA have similar results, and they perform worse than CCBBPSO-FID.
In Table V, the results show that CCBBPSO-FID always obtains the minimum cost results on all the tested instances, and only its mean result is lower than 45 000.Among the three well-performed CC algorithms, both CCDE-FID and CCBBPSO-AD are worse than CCBBPSO-FID.CSO, the algorithm for large-scale optimization, performs worse than CCBBPSO-FID and obtains better results than GPSO, BBPSO, CCDE-FID, and CCBBPSO-AD.Besides, SEO and RDA cannot obtain the feasible solutions in most test groups.
As shown in Table VI, although the data scale increases, the performance of CCBBPSO-FID does not deteriorate, and all results of it on Test III are still the best.Other versions of the CC algorithms, CCDE-FID and CCBBPSO-AD, perform better than other contestant algorithms.The results of CSO are better than those of classical PSOs (GPSO and BBPSO) and are beaten by the CC algorithms.In conclusion, CCBBPSO-FID shows an absolute advantage over the compared algorithms in different data scales.In addition, with the data scale increasing, CCBBPSO-FID still remains the best algorithm.Generally speaking, CSO, the popular algorithm for the large-scale optimization problems, also has promising performance on the LUSCND problem, and often does much better than GPSO, BBPSO, SEO, and RDA.However, the CC including CCDE-FID, CCBBPSO-AD, and CCBBPSO-FID, show a better performance on Test II and Test III.It indicates that the CC strategy is more effective to deal with the LUSCND problem and more suitable for large-scale optimization.The excellent performance of both CCDE-FID and CCBBPSO-FID further shows that the FID is a promising strategy for largescale optimization.Besides, CCDE-FID is outperformed by CCBBPSO-FID.The reason may be that some promising solutions are eliminated in the selection operator of DE, which the population diversity especially for the complicated COP.However, BBPSO reserves these solutions even if they are not improved.
For further illustration, box charts of three instances are shown in Fig. 10, and situations of other instances are similar.I_0_YL_DL_Norm, II_0_YL_DL_Norm, and III_0_YL_DL_ Norm are with the uncertain setting of YL, DL, and Norm, and are chosen from Test I, Test II, and Test III, respectively.As revealed by figures, CCBBPSO-FID not only performs much better than contestant algorithms but also is much more robust (with the volume).It indicates that CCBBPSO-FID is solve the LUSCND problem of smallscale or large-scale data.Moreover, the Wilcoxon rank-sum test at a 5% significance level is used for the statistical comparisons.The results of nine typical instances from three different data scales are shown in Table S.VI in the supplementary material.It can be seen that CCBBPSO-FID has a significant advantage over all the other compared algorithms on the nine instances.
2) Comparisons Regarding the Sensibility to Uncertainties: From the aspect of S_yield, based on the observation of From the aspect of C_demand, it can be observed that results are smaller when the demand level is higher (from DL to DH). High demand level means that there is more orders of suppliers.As a result, the unused capacity of suppliers will be less, and the corresponding penalty will decrease.Results are also different with different PDFs of demands (Norm, Log, and Tri).From the observation of Tables V and VI, it can be concluded that the results of Norm are smaller than those of Log and Tri, and the results of Log are the biggest.
To further analyze the influence of uncertainties on these algorithms, the mean maximum difference between results of different uncertainties is recorded, as shown in Table VII.For example, for the results of CCBBPSO-FID on Test I, the maximum difference among YL/YM/YH (6806 − 6701 = 105), DL/DM/DH (6760 − 6746 = 14), Norm/Log/Tri (6792 − 6716 = 76) and is 105.It can be seen that the difference of CCBBPSO-FID is smaller than that of other algorithms, which reveals that CCBBPSO-FID is less sensible to uncertain factors.The slight sensitivity is mainly due to that FID divides solutions into relatively independent parts which are solved by different subpopulations, helping to reduce the influence of uncertain factors on the whole.The slight sensibility also helps to enhance the robustness of CCBBPSO-FID, as shown in Fig. 10.

D. Analysis of Convergence Speed
To analyze the convergence speed of algorithms, the convergence curves are drawn for three instances, as shown in Fig. 11, and situations of other instances are similar.Since the maximum FEs of algorithms are different within the same execution time, the length of the line for each algorithm is different in the figures.From the observation of figures, it can be concluded that the convergence speeds of the CC algorithms (CCBBPSO-FID, CCBBPSO-AD, and CCDE-FID) are slower than those of other algorithms in the beginning.However, after some iterations, the evolution in non-CC algorithms (GPSO, BBPSO, CSO, SEO, and RDA) stagnates quickly.On the other hand, the CC algorithms surpass other algorithms afterward, except for CCBBPSO-AD on I_0_YL_DL_Norm, and CCBBPSO-FID obtains the best results in the end.The reason is that the CC algorithms solve the problem with four subpopulations and take four times FEs more than other algorithms to optimize the entire solutions.Therefore, they converge slowly at first but then converge quickly with the cooperative cooperation and the accumulated evolutionary information of all subpopulations.Moreover, FID decomposes the problem effectively based on different functions of different SCN members in the LUSCND problem, and it also helps CCBBPSO-FID and CCDE-FID to obtain the promising solutions.

E. Effectiveness of the CC Strategy
By comparing BBPSO with CCBBPSO-FID and CCBBPSO-AD in Tables V and VI, it can be seen that with the CC strategy, two CCBBPSO algorithms obtain better results, and the difference is obvious.CCBBPSO decomposes the problem into different parts, and these parts are solved separately by different subpopulations.These operations enhance the local search ability effectively.However, BBPSO only uses one population to deal with the problem, which limits the search ability of the algorithm.Further, the CC algorithms also perform superior to other competitor algorithms.

F. Effectiveness of FID
Random decomposition and average decomposition are designed to validate the effectiveness of FID.As mentioned above, CCBBPSO-RD that decomposes the problem aimlessly cannot obtain feasible solutions in all tests.By analyzing Tables IV-VI, it can be revealed that CCBBPSO-AD shows poorer performance than CCBBPSO-FID.The reason is that the average decomposition only divides the problem in the same size without considering the characteristics of the problem structure.Differently, FID decomposes the problem based on different functions of SCN members, and these functions are used to deal with different parts of the solution space.As a result, FID helps the algorithm to solve the problem more effectively.In addition, it improves the overall convergence speed of CCBBPSO, as shown in Section IV-D.The execution number and the success number of LRO ap are recorded during the evolution.The mean success rate of LRO ap in all algorithms is shown in Table IX.Remarkably, CCBBPSO-FID has the highest success rate among all the algorithms.It indicates that solutions obtained by CCBBPSO-FID are much more promising and closer to feasible ones.Furthermore, with the data scale bigger, it is more difficult for compared algorithms to repair solutions successfully, but the situation is contrary to CCBBPSO-FID.It can also be concluded that compared with the average decomposition of CCBBPSO-AD, FID of CCBBPSO-FID is more capable of finding feasible solutions.The highest mean success rate of all instances is up to 7.46%, and it proves that LRO ap is effective for the LUSCND problems.

G. Effectiveness of LRO ap
V. CONCLUSION SCND has important commercial values because it can help enterprises manage and plan the production line well and, then, cut cost and increase revenue.LUSCND is a more popular variant that extends the SCND model with large-scale and uncertain factors, but it is more difficult for existing approaches to solve.In this article, the efficient CCBBPSO-FID algorithm is proposed to handle this complex optimization problem.For the uncertain factors, MCM is used to calculate the fitness values in the uncertain environment.For the large-scale issue, the solution structure is first modified with the dimensionality reduction.Then, the CC strategy is combined with BBPSO to decompose the problem and to solve the decomposed parts with different subpopulations.Moreover, FID based on different functions of SCN members is devised

Fig. 3 .
Fig. 3. Example of the LUSCND problem. the managers of SCND, Fig. 3 shows an example of how to use the solutions obtained by an optimization algorithm to the LUSCND problem at a period.In this solution, there are one open supplier S 1,3 (the third location is chosen for S 1 ), two open warehouses W 1,3 and W 2,2 (the third location is chosen for W 1 , and the second location is chosen for W 2 ), and four customers (C 1 , C 2 , C 3 , and C 4 ).The capacities of S 1,3 , W 1,3 , and W 2,2 are shown at the top left in the figure, and the initial inventory quantities of W 1,3 and W 2,2 are shown in parentheses.In this period, it is assumed that the yield rate of S 1,3 is 0.8, and WS_order and WC_ship (the decision variables) are shown as solid lines with a number on the left.The example can be stated in four steps.First, customers send demands to warehouses.It is assumed that the total demand quantities of the four customers are 3, 2, 4, and 3, respectively.It is noted that there is no need to calculate the demand quantity of each customer to each warehouse in the FE, since the related cost (the penalty from customers) can be obtained by WC_ship in (6) and(8).Second, warehouses with the optimal locations are chosen to serve the customers (W_open, W 1,3 and W 2,2 in this example), and send orders to the supplier (WS_order).For example, W 1,3 sends 2 units of order to S 1,3 .Third, the supplier with the optimal location is chosen (S_open, S 1,3 in the example) and transports goods to warehouses (S_yield × WS_order), shown as the dotted lines with a number on the right.The yield rate of S 1,3 is 0.8, and it can deliver 2 × 0.8 = 1.6 units of goods to W 1,3 .Similarly, S 1,3 sends 3 × 0.8 = 2.4 units of goods to W 2,2 .After getting the supply, the quantities of W 1,3 and W 2,2 including supplies and inventory are 1.6 + 3 = 4.6 and 2.4 + 2 = 4.4, respectively.Finally, the warehouses that have been chosen in the second step deliver goods to customers (WC_ship).As shown in the parentheses, C 1 , C 2 , C 3 , and C 4 get 3, 2, 2.6, 1.4 units of goods, respectively.The demands of C 3 and C 4 are not satisfied.

Fig. 5 .
Fig. 5. Calculation of the fitness values.(a) Calculation of a feasible solution.(b) Calculation of an infeasible solution.

Fig. 7 .
Fig. 7. Two examples of LRO.(a) Solution is feasible after repair.(b) Solution is still infeasible after repair.
medium, and high level of uncertainties for the yield rate of suppliers, respectively.Moreover, these values obey the normal distribution.C_demand is generated in three different probability distribution function (PDF) like normal (Norm), log-normal (Log), and triangular (Tri), and each kind of distribution includes three levels according to the coefficient of variances (CV) (DL: CV = 0.05; DM: CV = 0.1; and DH: CV = 0.2), where DL, DM, and DH mean low, medium, and high level of uncertainties for the demand of customers, respectively.For the Norm and Log distributions, their mean values are set randomly within [capacity, 2 × capacity].For the Tri distribution, its lower limit a, upper limit b, and mode c are set to 0, 1/capacity, and 2/(3 × capacity), respectively, to ensure the height of the triangular 2/(b − a) is 2 × capacity.Totally, there are 3 × 3 × 3 = 27 uncertain factors for each instance.Data at three scales are generated with the same values of S_n and W_n (S_n = W_n = 5), and other parameters are shown in TableII.For the fair comparison, the execution time of each algorithm is equal for the same test set.Each scale includes five instances, such as I_0, I_1, I_2, I_3, and I_4 for Test I.There are 3 (scales) × 5 (instances) × 27 (uncertainty factors) = 405 uncertain instances in total.

TABLE I NOTATIONS
IN THE LUSCND PROBLEM

W_open, xbest.WS_order, xbest.WC_ship.
(variables solved in swarm 1 ) and xbest.That is, those variables not solved in this swarm are borrowed from the global best solution xbest.

TABLE II SCALE
CONFIGURATIONS AND EXECUTION TIME OF THE INSTANCE SETS

TABLE VII MAXIMUM
DIFFERENCE BETWEEN RESULTS OF DIFFERENT UNCERTAINTIESTables V and VI, it can be revealed that results mostly become smaller as the yield level gets higher (from YL to YH).YH includes a wider range [0.6, 0.9] with smaller values of the yield rate.It can be seen that if the yield rate decreases, the corresponding transport and inventory fees will decrease.Therefore, results of YH are smaller than those of YM and YL.

TABLE VIII RESULTS
OF ALGORITHMS WITH AND WITHOUT LROAP ON TEST I Table VIII shows the results of algorithms with and without LRO ap on Test I.It is obvious that algorithms without LRO ap cannot obtain feasible results (the data are in italic type).Moreover, results of all algorithms without LRO ap on Test II and Test III are approximately equal to 10 000 000 and not shown.It means that all algorithms without LRO ap cannot obtain feasible solutions in large-scale data.