A coevolutionary algorithm based on the auxiliary population for constrained large-scale multi-objective supply chain network

1.   Introduction

In the traditional manufacturing industry, enterprises always focus on their own single business, such as supplying raw materials, producing parts, assembling parts, or transporting. The isolating management pattern may be uncompetitive in the face of the volatile market, because enterprises cannot obtain complete information in time [1]. Afterwards, the pattern of supply chain management is proposed. Supply chain network (SCN) organizes enterprises from different stages together, and these enterprises work cooperatively and create a win-win situation [2,3]. Therefore, supply chain network (SCN) plays a significant role in the whole operation of the enterprise [4,5].

In the practical application of supply chain network, researchers can solve the SCN optimization problem through the mathematical modeling, the simulation, the data-driven, the game model, and other methods [6–8]. Evolutionary algorithms adopt a set of solutions (i.e., the population composed of multiple solution individuals) to solve problems, which have good global optimization ability and have become important methods to solve complex optimization problems in recent years [9,10]. For example, to solve the constrained optimization problems, evolutionary algorithms were combined with the random direction repair which was an effective constraint-handling method [11]. A constraint-objective cooperative coevolution framework was proposed to solve the large-scale constrained optimization problems, which allocated different computing resources to the decomposed sub-problems according to their contributions [12]. Therefore, evolutionary algorithms are often used to solve the complicated SCN optimization problems [13–15].

For example, Wang et al. [16] used genetic algorithm to choose suppliers and distribution centers and optimize the production quantity and the transportation volume in the iron and steel supply chain network, so as to minimize the operation cost and carbon emission of the network system. The problem is a multi-objective optimization problem of mixed integer programming with constraints. Sun et al. [17] used ant colony optimization algorithm to solve the scheduling optimization problem of mass customization supply chain and to minimize the production scheduling time by determining the optimal scheduling of cooperative suppliers. Zhang et al. [18] used the cooperative particle swarm optimization algorithm to solve the large-scale supply chain network design problem in the uncertain environment and to minimize the operation cost of the whole system by determining the optimal supplier and warehouse location scheme and the transportation logistics between network nodes. Akkad et al. [19] used a multi-objective heuristic approach to solve the collection and distribution problems of city logistics and to minimize both the fuel consumption and the emission of greenhouse gases. Saragih et al. [20] used a heuristic method to solve a location-inventory-routing problem in a three-echelon supply chain system with large scales. These problems have the characteristics of large scale, multiple objectives, or constraints which are the common difficulties of the SCN problems in the reality. However, the researchers did not solve the problem which has the three characteristics simultaneously.

In this paper, a constrained large-scale multi-objective supply chain network (CLMSCN) model is proposed, which is much closer to the reality and considers large scale, multiple objectives, and constraints in the same time. This model is to optimize the total operation cost and the customer satisfaction under the capacity constraints simultaneously. The total operation cost includes the production cost of the products, the inventory cost of the products, and the transportation cost among different SCN members. The proposed CLMSCN is a demand-driven model, which starts from the customers raising the demand to the SCN platform. According to the demand quantity of customers, the SCN members on the platform, including distributors and suppliers, begins to ordering, supplying, or distributing. This model is more competitive in the industry 4.0 era, since it starts from the market demand directly and helps to improve the operation efficiency with lower cost.

Besides, this paper proposes a coevolutionary algorithm based on the auxiliary population (CAAP), which uses two populations to solve the CLMSCN problem with complicated constraints. The first population is to solve the original problem with constraints, and a linear repair operator is designed to improve the feasibility of the infeasible solutions. The second population is to solve the problem without any constraints, which plays an auxiliary role to explore more solution space. With the cooperative work of the two populations, the CAAP algorithm can solve the CLMSCN problem in the experiment.

The rest of this paper is organized as follows. Section 2 describes the proposed CLMSCN problem in detail. Section 3 introduces all the components of the proposed CAAP algorithm. Afterwards, Section 4 discusses and analyzes the experimental results from different angles. Finally, Section 5 gives the conclusion of this paper.

2.   Problem description of CLMSCN

To optimize the total operation cost and the customer satisfaction, the CLMSCN model is proposed in this paper. The illustration of the CLMSCN problem is shown in Figure 1, which involves three kinds of members (including suppliers, distributors, and customers) and four stages (including customers sending demands to distributors, distributors sending orders to suppliers, suppliers supplying materials to distributors, and distributors transporting materials to customers). The whole process will last for several time periods. The numbers of suppliers, distributors, customers, and time periods are denoted as S, D, C and T, correspondingly. Before the description of the four stages of the problem, all variables involved in the CLMSCN problem model are described as follows:

Indices:

i: index of suppliers, i ∈ {1, …, S}

j: index of distributors, j ∈ {1, …, D}

k: index of customers, k ∈ {1, …, C}

t: index of time periods, t ∈ {1, …, T}

Parameters:

S: number of suppliers

S_cap_i: capacity of the supplier i

S_pdCost_i: production cost per unit capacity and time of the supplier i

D: number of distributors

D_cap_j: capacity of the distributor j

D_inv_{j, t}: inventory quantity of the distributor j in the period t

D_inv_{j, 0}: initial inventory quantity of the distributor j

D_invCost_j: inventory cost per unit capacity and time of the distributor j

C: number of customers

C_dem_k, _t: demand quantity of the customer k in the period t

SD_tpCost_{i, j}: transport cost per unit capacity and time from the supplier i to the distributor j

DC_tpCost_{j, k}: transport cost per unit capacity and time from the distributor j to the customer k

T: number of time periods

Decision variables:

DS_ord_{i, j, t}: order quantity of the supplier i from the distributor j in the period t

DC_tp_{j, k, t}: transport quantity of the distributor j to the customer k in the period t

The whole process of the CLMSCN problem in Figure 1 can be described as follows:

Stage 1) Demand: Each customer sends demands to the supply chain platform in each time period (C_dem_k, _t). It should be noted that the total demand quantity of the customer k is known, but the demand quantity of the customer to each distributor is unknown.

Stage 2) Order: If the previous inventory cannot satisfy the customer demands, distributors will send orders to the suppliers (DS_ord_{i, j, t}). In this stage, two capacity constraints should be satisfied. Firstly, in the period t, the total order quantity and the previous inventory of the distributor should not be larger than its capacity (Σ_i(DS_ord_{i, j, t}) + D_inv_{j, t-1} ≤ D_cap_j). Secondly, the total order quantity of a supplier should not be larger than its capacity (Σ_j(DS_ord_{i, j, t}) ≤ S_cap_i). Besides, the production cost of suppliers should be expended (Prod_cost = Σ_iΣ_jΣ_t(S_pdCost_i × DS_ord_{i, j, t})).

Stage 3) Supply: Suppliers transport materials to distributors. The transport cost from suppliers to distributors (Tp_cost1 = Σ_iΣ_jΣ_t(SD_tpCost_{i, j} × DS_ord_{i, j, t})) and the inventory cost of distributors (Inv_cost = Σ_jΣ_t(D_invCost_j × D_inv_{j, t})) should be expended. Besides, the inventory quantity of distributors should be updated (D_inv_{j, t} = D_inv_{j, t-1} + Σ_i(DS_ord_{i, j, t}) − Σ_k(DC_tp_{j, k, t}) = Σ_l(Σ_i(DS_ord_{i, j, t}) − Σ_k(DC_tp_{j, k, t})) + D_inv_{j, 0}, l ∈ {1, …, t}).

Stage 4) Transport: Distributors transport materials to customers (DC_tp_{j, k, t}). The transport cost from distributors to customers (Tp_cost2 = Σ_jΣ_kΣ_t(DC_tpCost_{j, k} × DC_tp_{j, k, t})) should be expended. When the demands of all customers are completely satisfied, the customer satisfaction is highest.

There are two assumptions in the problem model:

1)  The production quantity of a supplier is equal to the total order quantity.

2)  The production time of suppliers and the transport time from suppliers to distributors and from distributors to customers are ignored.

The two objectives of this problem is to minimize the operation cost (f1 = Prod_cost + Tp_cost1 + Inv_cost + Tp_cost2) and to maximize the customer satisfaction (f2 = Σ_kΣ_t(Σ_j(DC_tp_{j, k, t})/C_dem_k, _t)), which are mutually exclusive. The second objective can be equivalent to minimize f2 = Σ_kΣ_t(C_dem_k, _t)/ (Σ_jΣ_kΣ_t(DC_tp_{j, k, t}) + 1.0).

Based on the above descriptions, the final mathematical model of this CLMSCN problem can be described as follows:

Minimize

Subject to

The operation cost (f1) and the customer satisfaction (equivalent to f2) are calculated as Eqs (1) and (2), correspondingly. As for the domain of variables, the decision variables DS_ord_{i, j, t} and DC_tp_{j, k, t} and the inventory quantity of distributors (D_inv_{j, t}, which is updated in the third stage Supply) are all not smaller than zero, as shown in Eqs (3), (4), and (6), correspondingly. As for the capacity constraints in the second stage (Order), for the supplier i, its order quantity from all distributors (Σ_j(DS_ord_{i, j, t})) should not exceed its capacity (S_cap_i), as shown in Eq (5). For the distributor j, its order quantity to all suppliers (Σ_i(DS_ord_{i, j, t})) and previous inventory (D_inv_{j, t-1}) should not exceed its capacity (D_cap_j), as shown in Eq (7).

3.   CAAP algorithm

3.1.   Solution encoding

The solution in the proposed CAAP approach is encoded into a single dimensional array, as shown in Figure 2. Each solution has (dim = S×D×T+D×C×T) dimensions. The first (S×D×T) dimensions are consisted of DS_ord, and the last (D×C×T) dimensions are consisted of DC_tp.

3.2.   Initialization

Before the iteration, the two populations in the CAAP approach are randomly initialized. To generate more feasible solutions, reasonable lower and upper bounds should be determined. From the observation of Eqs (3) and (4), the lower bounds of DS_ord and DC_tp of the solution can be set as 0. For DS_ord of a solution, the decision variable DS_ord_{i, j, t} represents the order quantity of the distributor j to the supplier i in the period t, so it should not be larger than the capacity of the supplier i and the distributor j. Therefore, the upper bound of DS_ord_{i, j, t} is the smaller one between S_cap_i and D_cap_j (min(S_cap_i, D_cap_j)). For DC_tp of a solution, the decision variable DC_tp_{j, k, t} represents the transport quantity of the distributor j to the customer k in the period t, so it should not be larger than the capacity of the supplier i and the demand quantity of the customer k. Therefore, the upper bound of DC_tp_{j, k, t} is the smaller one between D_cap_j and C_dem_{k, t} (min(D_cap_j, C_dem_{k, t})). Finally, the decision variables are initialized randomly within the reasonable range.

3.3.   Linear repair operator

If an infeasible solution is generated in the population, the linear repair operator will be used to improve the feasibility of the solution. It can be seen that the constraints (Eqs (5)–(7)) of this problem are all the linear functions of decision variables (DS_ord and DC_tp). Therefore, to improve the feasibility of solutions, the decision variables can be linearly changed according to the constraints Eqs (5)–(7).

From the observation of the constraints Eqs (5)–(7), three conclusions can be obtained as follows:

1)  To improve the solutions which violate the constraint Eq (5), the decision variables in DS_ord should decrease.

2)  To improve the solutions which violate the constraint Eq (6), the decision variables in DS_ord should increase, or the decision variables in DC_tp should decrease.

3)  To improve the solutions which violate the constraint Eq (7), the decision variables in DS_ord should decrease, or the decision variables in DC_tp should increase.

Therefore, to maintain consistency as much as possible, three rules are formulated to repair infeasible solutions as follows:

1)  To repair the infeasible solutions violating the constraint Eq (5), the decision variables in DS_ord will decrease proportionally.

2)  To repair the infeasible solutions violating the constraint Eq (6), the decision variables in DC_tp will decrease proportionally.

3)  To repair the infeasible solutions violating the constraint Eq (7), the decision variables in DS_ord will decrease proportionally.

For example, for the infeasible solutions violating the constraint Eq (5), the new DS_ord_{i, j, t} (j ∈ {1, …, D}) need to approximately decrease by (Σ_j(DS_ord_{i, j, t}) − S_cap_i)/D proportionally, and the details are shown in Algorithm 1. For the infeasible solutions violating the constraint Eq (6), DC_tp_{j, k, t} (k ∈ {1, …, C}) need to approximately decrease by –D_inv_{j, t}/C proportionally, and the details are shown in Algorithm 2. For the infeasible solutions violating the constraint Eq (7), DS_ord_{i, j, t} (i ∈ {1, …, S}) need to approximately decrease by (Σ_i(DS_ord_{i, j, t}) + D_inv_{j, t-1} – D_cap_j)/D proportionally, and the details are shown in Algorithm 3.

It should be noted that the repair operator is time-consuming. Therefore, a self-adaptive repair probability of each solution is set. The initial repair probability of each solution is set as 10^-4 × dim. After an evaluation, if the solution is repaired, its repair probability will decrease by 10^-5 × dim; otherwise, the probability will increase by 10^-5 × dim.

3.4.   Complete CAAP algorithm

The CAAP algorithm applies the coevolutionary constrained multi-objective optimization framework [21], and uses two populations (pop1 and pop2) to solve the CLMSCN problem. pop1 is to solve the original problem with complex constraints, and pop2 is to solve the CLMSCN problem without constraints. Besides, the linear repair operator is designed to repair the infeasible solutions in pop1 as much as possible.

Algorithm 4 shows the pseudo code of the CAAP algorithm in details. The input parameter, Popsize, is the population size of pop1 and pop2, and the output parameters, PS and PF, are the Pareto solutions (the best solutions of a multi-objective problem) and their fitness values correspondingly. Firstly, pop1 and pop2 are initialized randomly, and the infeasible solutions in pop1 are repaired by the linear repair operator. Then, the iterative process is conducted, and the termination condition of this algorithm in line 4 is set as the maximal running time. In the iterative process, Parent1 and Parent2 are selected from Parent2 correspondingly by the roulette wheel selection. Then, in lines 7 and 8, offsprings are generated by NSGA_II with the simulated binary crossover and polynomial mutation [22–24]. Then, pop1 and pop2 are updated and evaluated in lines 9 to 13. In lines 12 and 13, the fast-non-dominated-sort in [22] and the farthest-candidate approach in [25] are used to select the non-dominated solutions. Finally, PS and PF are updated.

4.   Experimental verification and comparisons

4.1.   Experimental settings

To validate the performance of the CAAP algorithm, three different scales of instance sets are randomly generated. Each scale has five instances, such as I_1 to I_5 of the small scale, II_1 to II_5 of the middle scale, and III_1 to III_5 of the large scale, and the configurations and the maximum execution time of the instance sets are shown in Table 1. Besides, MPCMO_BBPSO (multiple populations for constrained multi-objective optimization with bare-bones particle swarm optimization), PBBPSO (bare-bones particle swarm optimization with the Pareto optimality), NSGA_II (nondominated sorting genetic algorithm II, [22]), MOEA/D (multi-objective evolutionary algorithm based on decomposition, [26]), and NSLS (nondominated sorting and local search, [25]) are used to compare with the CAAP algorithm for the CLMSCN problem. MPCMO_BBPSO, which applies the multiple populations for multiple objectives framework [27], takes the constraints as a new objective and uses the bare-bones particle swarm optimization to solve the CLMSCN problem. PBBPSO uses bare-bones particle swarm optimization which is competitive in solving the supply chain optimization problem. NSGA_II, MOEA/D, and NSLS are the competitive algorithms for the multi-objective optimization problems. For the fair comparison, the population size of all algorithms is set 50, and the linear repair operator proposed in this paper is also used in the all compared algorithms to repair the infeasible solutions.

4.2.   Compared results

To evaluate the performance of the CAAP algorithm and other contestant algorithms from diverse angles, two metrics are used, HV (hypervolume) and C(A, B) (A and B are two algorithms) [28]. The HV value represents the hypervolume surrounded by the reference point and the Pareto front. For the minimization problem solved in this paper, if the HV value of an algorithm is higher, it means the algorithm is better than other algorithms. The C(A, B) value means the ratio of the solutions which are obtained by the algorithm B and are dominated by the solutions obtained by the algorithm A. Therefore, if the C(A, B) value is larger than the C(B, A) value, it means that the algorithm A can obtain better solutions than the algorithm B.

Table 2 shows the HV values of all algorithms on different instances. The data in bold represent the best result among all algorithms. It can be observed from Table 2 that for the instances in small scale, CAAP can obtain the best results on 3 instances, followed by MPCMO_BBPSO and MOEA/D. For the instances in middle and large scale, CAAP can obtain the best results on all instances, followed by MPCMO_BBPSO and NSGA_II. With the problem scale increasing, the advantage of the performance of CAAP does not deteriorate.

Table 3 and Table 4 show the C(A, B) values of all compared algorithms with CAAP on different instances. From the observation of Table 3, CAAP can obtain 2, 5 and 4 better results than MPCMO_BBPSO, PBBPSO, and NSGA_II on small instances, correspondingly. For the middle and large instances, CAAP obtains the best results on all instances compared with MPCMO_BBPSO, PBBSO, and NSGA_II. From the observation of Table 4, CAAP outperforms MOEA/D and NSLS on all different instances.

For further illustration, Figure 3 depicts the Pareto fronts of all algorithms on three instances from different problem scales. For the minimization problem solved in this paper, if the Pareto front of an algorithm is the closest to the coordinate axes, it means that the algorithm can get the best results among these algorithms. Besides, if the Pareto front of an algorithm is longest, it represents that the algorithm can obtain most non-dominated solutions that spread widely in the solution space. It can be observed from Figure 3 that the Pareto fronts of CAAP on the three instances are almost all the closest and longest among all algorithms, which demonstrates that CAAP can obtain the best results on different instances for the CLMSCN problem.

From the above observations, CAAP can outperform the compared algorithms on different instances. The reason may be that CAAP uses two populations to cooperatively solve the complicated constrained multi-objective problems. One population is to solve the original problem with complex constraints, and the linear repair operator with the self-adaptive repair probability is used to improve the feasibility of solutions. The other population is to solve the simplified problem without constraints, which helps to explore more solution space quickly. The cooperative mechanism of the two populations is effective to solve the complex constrained problem.

4.3.   Effectiveness of the linear repair operator

To validate the effectiveness of the linear repair operator, the repair probability of all algorithms on different instances are recorded in Table 5. It can be observed that the repair probability of all algorithms is all much larger than zero, and with the problem scale increasing, the repair probability increases. The results illustrate that the linear repair operator can help algorithms improve the feasibility of solutions. Besides, the C(A, B) of CAAP and CAAP_noRep (CAAP without the linear repair operator) on different instances is shown in Table 6. It should be noted that CAAP_noRep uses the superiority of feasible solutions [29] as the constraint handling method. CAAP can get much better results than CAAP_noRep, especially on middle and large instances. It can also be concluded that the linear repair operator is effective to improve the solution quality of algorithms for solving the CLMSCN problem.

5.   Conclusions

This paper designs the CLMSCN model with complex constraints and multiple objectives which both considers the minimization of the total operation cost and the customer satisfaction as the optimization objectives. Besides, to solve the complex problem effectively, the CAAP algorithm is proposed in this paper, which applies two populations to solve the problem. The first population is to solve the original CLMSCN problem with complex constraints, and the linear repair operator will be used to improve the feasibility of solutions. The second population is to solve the simplified problem without constraints, which helps to quickly expend the search space of solutions. The experimental results show that the cooperative mechanism of CAAP helps to obtain much better solutions than the compared algorithms on the CLMSCN problem with complex constraints and multiple objectives, especially on the large-scale instances. In the future work, we will design more relatively general constraint-handling methods for the complex problems, and apply the cooperative mechanism of CAAP to solve other supply chain network optimization with complex constraints and multiple objectives, such as the minimization of the fuel consumption and the emission of greenhouse gases.

Acknowledgments

The work described in this paper was substantially supported in part by the National Key R & D Program of China under Grant 2017YFC1601000 and 2017YFC1601800, the National Natural Science Foundation of China under Grant No. 62106088, No. 62172192, and No. 62073155, in part by the High level personnel project of Jiangsu Province (JSSCBS20210852), and in part by the Fundamental Research Funds for the Central Universities, Jiangnan University and Jilin University.

Conflict of interest

The authors declare there is no conflicts of interest.