MEA-CNDP: A Membrane Evolutionary Algorithm for Solving Biobjective Critical Node Detection Problem

The critical node detection problem (CNDP) refers to the identification of one or more nodes that have a significant impact on the entire complex network according to the importance of each node in a complex network. Most methods consider the CNDP as a single-objective optimization problem, which requires more prior knowledge to a certain extent. This paper proposes a membrane evolution algorithm MEA-CNDP to solve biobjective CNDP. MEA-CNDP includes a population initialization strategy based on the evaluation of decision variables, a strategy to transform the main objective, a strategy to update the membrane inherited pool, and four membrane evolutionary operators. The numerical experiments on 16 benchmark problems with random and logarithmic weights show that MEA-CNDP outperforms other algorithms in most cases. In particular, MEA-CNDP has unique advantages in dealing with large-scale sparse bi-CNDP.


Introduction
e nature-inspired optimization strategy has shown very interesting results in many fields, such as teach optimization techniques and so on [1]. As a new branch of natural computing, membrane computing has attracted continuous attention of researchers, it is a computing model abstracted from biological cells, organs, and tissues. Based on the cell membrane, the membrane computing establishes a calculation model based on the function and characteristics of the cell membrane. Such a calculator is also called the P system. Membrane computing has been widely concerned since it was proposed, and it has shown excellent results in many applications. Most P systems have been proved to have the same computing power as Turing machines [2,3]. In theory, many NP-hard problems can also be solved in polynomial time using the P system, such as the all-set problem [4] and Hamilton loop problem [5]. Since membrane computing was put forward, it has received extensive attention and has been applied in many fields, such as image processing [6], numerical optimization [7], and social network analysis [8].
Inspired by membrane computing, Nishida combined membrane nested structure with other evolutionary algorithms in 2004 and proposed the membrane algorithm framework (MA) for the first time [9]. Since then, researchers call the algorithms based on this framework membrane algorithm, which provides a nested structure compatible with other algorithms, making MA well combined with various algorithms to solve different problems. For example, a membrane clustering algorithm for automatic clustering is proposed in [10] based on MA, which can easily find the number of clustering under the control of evolutionary communication mechanism. Combining particle swarm optimization (PSO) and MA, a membrane heuristic algorithm based on particle swarm optimization (mMPSO) is proposed to solve the path planning problem of multiobjective robots in dynamic obstacles and dangerous environments [11]. Moreover, in [12], MA is combined with differential evolution algorithm (DE), in [13], and MA is combined with quantum evolution algorithm (QEA).
Based on MA, membrane evolutionary algorithm (MEA) is introduced, MEA is an evolutionary algorithm framework inspired by the behavior characteristics of biological cells. It simulates various life activities of cells, realizes information exchange between cells by relying on membrane, and thus realizes individual evolution. In [14], a membrane evolutionary algorithm for solving minimum vertex cover problem (MEAMVC) is introduced, and the experimental results show that the MEAMVC has good performance in solving the minimum vertex coverage problem. In [15,16], MEA are introduced to solving the maximum clique problem and travelling salesman problem, respectively.
A crucial research direction in engineering design is to consider optimization problems. With the continuous increase of the scale of optimization problems and the increase of dynamic and random factors, traditional methods reflect the complexity and instability of calculation. erefore, researchers are looking for various new methods to solve optimization problems. In recent years, some excellent methods have emerged and have been applied to various fields. For example, Zapata et al. [17] propose that a novel hybrid swarm algorithm combining strengths of self-assembly and the particle swarm optimization greatly improves the convergence speed. Reference [18] presents a novel application of the metaheuristic slime mould algorithm (SMA) to the optimal tuning of interval type-2 fuzzy controllers. Population-based evolutionary algorithm has shown gratifying results in solving optimization problems, and a new estimation of distribution algorithm (EDA) is introduced in [19]. It maintains a trade-off between run time and the number of evaluation points. Reference [20] combines the island model and cuckoo search (CS) algorithm and thus presents the islandbased CS with polynomial mutation (iCSPM), which greatly improves the search ability and maintains the population diversity.
Complex networks can be seen everywhere in both scientific research and daily life. As an important subject of optimization problems, the critical node detection problem (CNDP) aims to identify the critical nodes in complex networks, which is also a common means to study graph problems. Similar problems include the critical node problem (CNP), which aims to minimize the connectivity of the residual graph under the constraint of the maximum number of nodes that can be deleted [21]. Cardinalityconstrained critical node detection problem (CC-CNP) is devoted to minimizing the number of deleted nodes given the maximum connected component [22].
At present, the methods to solve CNDP mainly include accurate algorithm and heuristic algorithm. For the accuracy strategy, most algorithms consider CNDP as an integer programming problem [23][24][25], which has significant limitations. In contrast, the heuristic algorithm has achieved good results in solving CNDP. A hybrid heuristic algorithm is proposed with the combination of greedy algorithm and local search algorithm and shows good performance in [26]. In [27], incremental learning based on population is combined with a simulated annealing optimization algorithm based on combinatorial rank free problem representation to solve CNDP, which saves space and reduces the need for individual repair mechanism simultaneously. Most of the above algorithms treat CNDP as a singleobjective optimization problem. Among them, the maximum number of nodes that can be deleted and the maximum number of connected components are often required. However, in practical problems, prior knowledge is often not easy to obtain. Moreover, considering the network connectivity and the cost of deleting nodes are often two conflicting objectives, which should be viewed simultaneously. erefore, establishing a biobjective optimization model for CNDP is considerable. In [28], CNDP is considered a biobjective optimization problem for the first time, and the Pareto fronts of some problems are given. In [29], a specific biobjective optimization model of CNDP is proposed, and the ant colony algorithm is used to solve it. In [30], a biobjective optimization model for CNDP is established, and the two objective functions are the number of connected components and their cardinal variance, respectively, but this is actually a generalization of pairwise connectivity. is paper proposes a novel evolutionary algorithm based on MA for solving the biobjective critical node detection problem (bi-CNDP) model in [31], called MEA-CNDP for short. Based on membrane division, differentiation, death, and other life activities, MEA-CNDP designs and implements four evolutionary operators: division, fusion, cytolysis, and selection operator. Aiming at bi-CNDP, a new population initialization strategy is proposed in MEA-CNDP, and the main objective transforming and membrane-to-subproblem matching strategy are adopted to improve the efficiency of MEA-CNDP. e main contributions of this paper include the following: (i) Taking into account the characteristics of bi-CNDP, a new population initialization strategy based on the evaluation of decision variables is proposed. By calculating the number of nondominated front and the cost of nodes, the importance of each decision variable can be levelled out, thereby generating a high-quality initial population. (ii) Design and implement four evolutionary operators, including division, fusion, cytolysis, and selection operator. By the communication between membranes, those evolutionary operators can generate better individual membrane, eliminate the membrane with poor quality, and update the external archive membrane. (iii) A membrane evolutionary algorithm MEA-CNDP is proposed to solve the bi-CNDP. e experimental results on four different types of instances show that the MEA-CNDP algorithm has unique advantages in dealing with large-scale sparse bi-CNDP. e rest of this paper is organized as follows. Section 2 introduces the basic knowledge and related definition, including the concept of MA and MEA, CNDP, and bi-CNDP. Moreover, some algorithms solving bi-CNDP are introduced briefly. In Section 3, the framework of MEA-CNDP is proposed, and the details of related strategies for improving the efficiency are introduced; at last, the implementation of four types of evolutionary operators is detailed. In Section 4, 2 Computational Intelligence and Neuroscience based on four types of bi-CNDP instances, the comparative experiments are designed, and the analysis of the experimental results is given as well, which verified the performance of MEA-CNDP. Finally, the conclusion and future work are drawn in Section 5.

Preliminaries
is section provides the related research basis of this paper. Firstly, the related content of membrane computing is introduced, and the related knowledge of MEA is submitted. en, the related definition and formula of CNDP are presented briefly, thus providing the multiobjective optimization problem (MOP) model of CNDP, that is, the biobjective critical node detection problem (denoted as bi-CNDP). Finally, some related algorithms solving bi-CNDP are introduced briefly, which are used in the subsequent comparative experiments.

Membrane Algorithm and Membrane Evolutionary
Algorithm. As a biological computing model abstracted from the level of biological cells, membrane computing is committed to studying the function and characteristics of the cell membrane in organisms and applying them to optimization problems; such a calculator is also called a P system [32]. P system is mainly composed of its structure, material object, and evolutionary rules. e main body of the P system is to construct the membrane structure corresponding to specific problems by simulating the basic structure of the cell membrane. Evolutionary rules are the basis of material communication within and between membranes. Different rules are implemented to transform different information.
MA is an algorithm framework derived from membrane computing. MA provides a nested structure compatible with various other algorithms, enabling various excellent algorithms to solve the final problem in the framework of membrane computing, aiming at specific subproblems, giving full play to their respective advantages. e nested structure of MA is shown in Figure 1. e main structural feature of MA is a layer of nested membranes. e region between each two nested membranes contains subalgorithms. In different regions, the subalgorithms can be the same or different. ese subalgorithms usually solve different stages of the problem or solve different parts of the problem.
MEA is an evolutionary algorithm based on cell division, fusion, death, and other life activities, according to the characteristics of the cell membrane to achieve information exchange between cells. According to the life activities of cells and the information exchange mechanism between cells, MEA designs the operators of division, fusion, cytolysis, and selection and uses the membrane parallel optimization and heuristic information to realize individual evolution and further solve the problem. e algorithm framework of the membrane evolutionary algorithm is shown in Algorithm 1. e framework of MEA includes four parts: population and membrane structure initialization, membrane evolution, membrane repair, and final solution output. e primary evolutionary process of MEA is the four operators of division, fusion, cytolysis, and selection. e membrane population is constantly operating these four operators to make the whole population evolve in a better direction.
(i) Division operator: in MEA, one individual membrane divides into two submembranes, the substance in the original individual membrane is distributed into the two submembranes, and the original individual membrane is removed. e process of division operator is shown in Figure 2(a). (ii) Fusion operator: in MEA, the fusion operator merges two individual membranes into one membrane, and all the substances in the original two membranes enter into the newly formed membrane, while the original two membranes and their substances are removed. e process of the fusion operator is shown in Figure 2(b). (iii) Cytolysis operator: in MEA, the individual membrane with poor performance or without specific conditions will be dissolved, and the dissolved substances in the membrane can be deleted or released into the environmental membrane. (iv) Selection operator: in MEA, the membrane with higher fitness or better individuals should be retained, and the selection operator is used to update the excellent individuals in the membrane population evolution.
In MEA, each individual membrane is regarded as a candidate solution, and information exchange and synchronization are realized among multiple individual membranes; thus, the whole membrane population evolves in a better direction. Precisely because of such characteristics, MEA possesses remarkable parallelism and scalability, flexible coding, and excellent search capability. Computational Intelligence and Neuroscience 3

CNDP and
includes the promotion or mitigation of the entire complex network.
Let G � (V, E) be an undirected weighted graph, the number of nodes is |V| � n, and the number of edges is |E| � m. e goal of CNDP is to determine a subset R⊆V, so that the residual graph G(V\R) has the minimum number pairwise connectivity (i.e., the number of a pair of nodes connected by a path in a graph), where |R| ≤ K, K is the maximum number of the deleted nodes allowed. If the nodes of CNDP are weighted, that is, the weighted critical node detection problem (denoted as wCNDP) transforms the constraint to ω R ≤ ω, where the ω R is the sum of the weights of the deleted nodes and ω is the maximum weights of the deleted nodes allowed. More detailed about wCNDP can be seen in [30].
To solve the CNDP as a single-objective optimization problem, it is often required to give the maximum number of deleted points or determine the maximum number of connected components. However, in practical issues, the prior knowledge is often not easy to obtain, and considering the network connectivity and the cost of deleted nodes are often two conflicting objectives, which should be viewed simultaneously. erefore, it is of great value to establish a biobjective optimization model for CNDP. e model of bi-CNDP is introduced in [31], which considers the pairwise connectivity and the cost of deleting nodes simultaneously. Let G � (V, E, C) be a weighted undirected graph, V � v 1 , v 2 , . . . , v n is the set of nodes of G, E � (v i , v j )|i, j ∈ 1, 2, . . . , n { } is the set of edges, and the cost of deleting each node is given in C � c 1 , c 2 , . . . , c n . e first objective function of bi-CNDP is about pairwise connectivity, which can be calculated by where x ij � 1 if v i and v j are connected; otherwise, x ij � 0; it can be further normalized as It can be seen that the value range of nPWC(G) is [0, 1], obviously, nPWC(G) can be considered the measure of connectivity of the graph because nPWC(G) is proportional to the number of nodes connected through at least one path. e closer the nPWC(G) value is to 1, the closer the graph is to the connected graph.
After deleting nodes, except for minimizing the connectivity of the network, minimizing the cost of deleting these nodes is considerable. erefore, the second objective function of bi-CNDP is used to measure the cost of deleting the selected nodes R⊆V, which can be calculated by Input: problem description, related parameters Output: the final solution (1) Membrane structure initialization (2) Population initialization (3) while termination criterion not fulfilled do (4) (1) Division operator: a membrane is divided into two submembranes, and the substances in the original membrane are distributed in the two submembranes (5) (2) Fusion operator: two membranes fuse to form a membrane, and the substances in the former two membranes belong to the same membrane (6) (3) Cytolysis operator: dissolve several membranes from the population, and the substance in the dissolved membrane are dissolved as well (7) (4) Selection operator: select several membranes from the current membrane population to enter the next generation   Computational Intelligence and Neuroscience where c i is the weight associated with v i , v i ∈G c i is drawn to normalizing, which makes the value range in [0, 1] as well. erefore, the model of bi-CNDP can be described as minimize(nPWC(G(V\R)), nCost(R)), s.t.
where y � (y 1 , y 2 , . . . , y n ) is the decision variables, y i is the binary variable, y i � 1 if the node v i is deleted, and y i � 0 otherwise. x ij shows the connectivity of v i and v j , x ij � 1 if v i and v j belong to the same connected component, and x ij � 0 otherwise.
y � (y 1 , y 2 , . . . , y n ) is the decision variables to be optimize, y i values in {0, 1} and changes continuously with the implementation of the programme, until it evolves into a good enough solution. e first constraint ensures that the residual graph after deleting some nodes includes several connected components. e second constraint ensures that, for three nodes in the same connected component, any two nodes between them are also in the same connected component. e third constraint specifies the range of values for x ij and y i .
It can be seen that bi-CNDP defined above is NP-hard on the general graph. Although the Pareto front of the above model cannot be obtained for the time being for the general graph, for the CNDP whose graph structure is a tree, the algorithm with polynomial time can get the approximate Pareto front of the above model. erefore, a heuristic algorithm is needed for general graphs to find the approximate Pareto front of the above bi-CNDP.

Related Algorithms
e strategy based on decomposition is a practical approach and widely used method in multiobjective optimization evolutionary algorithm, such as the multiobjective evolutionary algorithm based on decomposition (MOEA/D) is proposed in [33]. MOEA/D decomposes the MOP into several scalar single-objective optimization subproblems and then optimizes the subproblems simultaneously. In MOEA/D, each subproblem is optimized only according to the information of its adjacent subproblems, which significantly reduces the computational complexity of each generation.
In MOEA/D, it is crucial to decompose the MOP into multiple scalar subproblems, and the common aggregation functions mainly include the weighted sum approach, Chebychev approach, and boundary intersection approach [34]. MOEA/D firstly introduces the decomposition method into MOP, making the idea of decomposition into the evolutionary algorithm. In addition, by decomposing the MOP into multiple scalar optimization subproblems for optimization, MOEA/D significantly reduces the difficulty of maintaining population diversity and fitness allocation.

DMOEA-εC.
Based on MOEA/D, the decompositionbased multiobjective evolutionary algorithm with the ε-constraint framework (DMOEA-εC) is proposed in [35]. DMOEA-εC uses the ε-constraint method to MOPs, it selects one objective function as the main objective, transforms other objectives into constraints, allocates upper bound vectors for each subproblem, and then optimizes multiple subproblems simultaneously for MOP; the ε-constraint formula can be defined as is the upper bound vector, and the index of the main objective is selected from all the objective functions randomly or assigned by the decision-maker. ρ is a tiny positive number, T are the ideal optimal node and ideal worst node, respectively.
DMOEA-εC proposes the strategy of main objective alternation, which exchanges the main objective function through a specific strategy or periodically and matches the solution with the subproblem again after the main objective function alternation. After a new feasible solution is generated by evolution, DMOEA-εC fits the feasible solution with the subproblem again so that the feasible solution newly generated can be utilized to the maximum extent.

I-MOEA/D and I-DMOEA-εC.
MOEA/D and DMOEA-εC perform well on continuous or discrete benchmark problems, but for bi-CNDP, a more special strategy is needed to make them adopt the characteristics of bi-CNDP. us, the new mating pool and replacement pool are designed in [31].
In [31], there are four types of strategies for mating pool, including selecting from the neighbors, the whole population, the neighbors and the population, and the external archive population.
ere are two strategies for the replacement pool, including the local replacement and the global replacement. e experimental results have determined the two algorithms with the best performance, namely, I-MOEA/D and I-DMOEA-εC.
It is worth noting that so far there are few special algorithms to solve bi-CNDP. As the variant of algorithm MOEA/D, I-MOEA/D, and I-DMOEA-εC has shown good results in solving bi-CNDP, inspired by this, this paper aims to explore new better algorithm to solve bi-CNDP; more details are introduced in Section 3.

The Membrane Evolutionary Algorithm Solving Bi-CNDP
In view of the bi-CNDP introduced above, this section designs and implements the membrane evolutionary algorithm MEA-CNDP solving bi-CNDP. Representation of candidate solutions and related definitions are presented at first. en, according to the characteristics of bi-CNDP, the membrane evolutionary algorithm framework of MEA-CNDP is designed. Next, some specific strategies are introduced, including population initialization strategy, membrane inherited pool updating, and main objective transforming strategy. Finally, the membrane evolutionary operators, including division, fusion, cytolysis, and selection, are introduced in detail.

Representation and Definition of Candidate Solutions.
e initial membrane structure of the MEA-CNDP algorithm is shown in Figure 3. Membrane 0 is the environment membrane (also known as skin membrane) and provides the foundation of the whole membrane population evolving. In the environment membrane, there are the membrane population (membrane 1, 2, . . ., N) for evolution, and the inherited pool (recorded as IP) for each iteration, the external archive membrane (registered as EM) for recording the nondominated solutions. Moreover, all decision variables and upper bound vectors keep only one copy of them in the environment membrane.
In the process of population evolving, when the individual membrane (membrane i) needs to fuse decision variables or upper bound vectors, the decision variables or upper bound vectors are copied into the individual membrane. When an individual membrane is divided, the decision variables and upper bound vectors in it will be dissolved, which means the number of decision variables and upper bound vector in the environment membrane will remain unchanged.
Each evolution only takes place in the membrane population, that is, the division, fusion, cytolysis, and selection operators are only performed in membrane 1, 2, . . ., N. Membrane IP is the current inherited pool, which stores the individual membranes participating in the following evolutionary process. Moreover, the membrane EM is used to record the nondominated solutions in the current membrane population. When a new membrane is generated, the membranes in EM will be updated synchronously. After the evolutionary process is over, the membranes in the EM will output the final solution.

Framework of MEA-CNDP.
is section introduces the framework of the proposed MEA-CNDP algorithm. MEA-CNDP performs division, fusion, cytolysis, and selection operators in each evolutionary process and maintains the convergence of the solution by saving and updating the optimal solution. Before each evolutionary process, MEA-CNDP adopts specific strategies to optimize its performance. In Algorithm 2, MEA-CNDP firstly generates N uniformly distributed upper bound vector ε 1 , ε 2 , . . . , ε N by dividing each objective axis with equal spacing Δ. Figure 4 shows the process of generating upper bound vectors through a 3-objective optimization problem model (m � 3), where the equal spacing Δ � 1/4 (q � 5), 25 upper bound vectors are generated by this way.
MEA-CNDP designs a population initialization strategy; lines 3 and 4 calculate the scores of the decision variables and generate the initial population according to the scores of decision variables, respectively. e initial population generated by the strategy has good characteristics, accelerates the convergence of the population, and maintains diversity. Algorithms 3 and 4 show the details.
Input: P (current membrane population), NoF (the mark of main objective function) Output: P (membrane population) (1) S ← ∅; (2) for each m in P do (3) Remove the upper bound vector of m and add it to S; (4) while S is not empty do (5) ε l ← Randomly select a upper bound vector from S; (6) for each m i in P do (7) y NoF+1 ← Calculate the value of membrane m i to the (NoF + 1)-th objective function; (2) for l � 1 to N do Remove the decision variables of membrane m i ; (19) Duplicate Computational Intelligence and Neuroscience population, respectively. Z * and Z nad are updated with the evolutionary process.
e NoF and Cur Gen in lines 9 and 10 mark the main objective and iteration counter, respectively. NoF ∈ {0, 1}, when NoF � 0, the first objective function is selected as the main objective; otherwise, the second objective function is selected. Cur Gen is used to record the number of iterations, and the algorithm ended when Cur Gen meets the maximum iteration number Max Gen .
Lines 12-22 are the evolutionary process of each generation in MEA-CNDP. In IP updating (Line 13), the individual membranes in the original IP are removed completely, and Algorithm 5 selects the new individual membranes into IP. In the main objective transforming (Lines 15-16), the upper bound vectors and individual membranes are rematched according to Algorithm 6. e evolution of IP includes membrane division (Algorithm 7) and updated Z * , membrane fusion (Algorithm 8), membrane cytolysis and selection (Algorithm 9), and updated Z nad . Figures 5-8 illustrate the procedure of MEA-CNDP. Figure 5 shows the membrane structure after one iteration, and the value of Cur Gen will be determined before the next iteration. If the Cur Gen can be divided by DRA interval , Algorithm 5 will be executed, and if the Cur Gen can be divided by INm, Algorithm 2 will be executed. Here, we assume that Algorithms 5 and 6 have been executed. Figure 6 shows the process and results of executing the division operator. In one evolution, according to Algorithm 7, membranes m j and m 3 are selected to perform the division operation, and membrane o is the generated offspring membrane. Figure 7 shows the process and results of performing the fusion operator.
According to Algorithm 8, the upper bound vector of membrane o and the substances in the neighbor membrane of membrane o are updated. Figure 8 shows the process and results of performing the cytolysis and selection operator. According to Algorithm 9, the membrane o interacts with the membrane in the external archive membrane EM, dissolves the bad membrane, and selects the retained membrane.
MEA designs evolutionary strategy according to the information exchange between membranes. According to the life activities of membrane and the information exchange mechanism between membranes, MEA generates its own unique division, fusion, dissolution, and selection operators. MEA uses membrane heuristic information and parallel optimization to realize individual evolution and solve problems. In particular, MEA-CNDP takes into account the characteristics of the bi-CNDP and introduces new strategies such as population initialization and new evolutionary operators, thereby establishing a new framework to solve bi-CNDP.

e Initialization Strategy.
e pro and cons of the initial population significantly affect the bi-CNDP, especially for the models similar to large-scale MOPs. us, MEA-CNDP introduces a population initialization strategy, which includes two parts: evaluating the scores of decision variables and generating the initial population. rough this strategy, not only can the quality of the initial population be improved but it can also provide evidence for deciding whether to retain or eliminate evolutionary genes.  (6) Remove membrane m in EM; (7) if m dominate o then (8) isDominated ← true; (9) if o_F1 � � m_F1 && o_F2 � � m_F2 then (10) isDominated ← true; (11) if isDominated is true then (12) Add membrane o to EM; (13) return EM; ALGORITHM 9: Cytolysis and Selection. e process of calculating the scores of decision variables can be seen in Algorithm 3, the coding method is used in [37], and some improvement is executed according to the characteristic of bi-CNDP. Firstly, since the decision variable of bi-CNDP is a binary value (if the node is deleted, the value of the corresponding decision variable will be set to 1; otherwise, it is 0), Algorithm 3 simplifies the coding process, adopts binary coding directly, and omits the real vector coding of the solution. Next, in bi-CNDP, the degree of a node is taken into account. Nodes with degree 1 are not included in the score calculating of the decision variables (in fact, their scores have been assigned to Max_Double). Deleting nodes with degree 1 in the graph has little effect on the pairwise connectivity of the graph. Especially in the large-  Figure 5: e membrane structure before the beginning of a certain evolution.  x 1 , x 2 , x 3 , x 7 , x 8 , x 9 ε p x 2 , x 3 , x 4 , x 6 , x 7 , x 8 ε q ε 1 , ε 2 , . . ., ε k-1 , ε k , ε k-1 , . . ., ε N x 2 , x 3 , x 5 , x 7 , x 8 , x 9 ε 4 x 1 , x 2 , x 3 , x 7 , x 8 , x 9 ε p x 2 , x 3 , x 5 , x 7 , x 8 , x 9 ε q x 2 , x 3 , x 5 , x 7 , x 8 , x 9 ε k Computational Intelligence and Neuroscience scale graph, the effect is almost negligible. On the contrary, it will increase the impact on "minimize deleting cost." Finally, because the nodes in bi-CNDP have weights related to the cost of deleting themselves, the final calculation results will be determined by the nondominated front number and the cost of deleting the corresponding node together, which can more accurately assess the importance of each decision variable.
Algorithm 4 shows the process of generating the initial population. Each individual membrane is generated by adding decision variables to the empty membrane, and the selected decision variables are determined by comparing their scores. e size of each membrane is set as rand() × D, which is seen as a random parameter. e decision variables with smaller score values are more likely to be selected because the smaller score values often mean better nondominated solutions.

Membrane Inherited Pool Updating and Main Objective
Transforming. Different individuals in the membrane population have different upper bound vectors and subproblems, which means that the calculation difficulty is also different. erefore, it is reasonable to assign different amounts of computational effort to different individuals according to the utility value in each generation of evolution. Algorithm 5 updates the IP through fixed iteration interval, and in each evolutionary process, only the individual membranes in IP participate in the evolution, which ensures the efficiency of the MEA-CNDP algorithm.
Algorithm 5 calculates the utility value of each individual membrane in Lines 6-9. Every time the IP is updated, all individual membranes in the original IP are removed. Line 14 selects several individual membranes to enter the IP through binary tournament selection, which means that the individual membrane with better utility value will be more likely to enter the IP. is process is repeated (N/5 − m − 1) times to select adequate individual membranes, and these individual membranes in the IP participate in the evolutionary process until the IP is updated next time. e details of parameter selection can be seen in [35].
Since transforming the main objective is adopted, the main objective function will convert with a fixed iteration interval. Among the current membrane population, the individual membranes that perform well to the original main objective function may not perform well to the current main objective function; thus, it is necessary to rematch the membrane and the subproblem.
Algorithm 6 matches the individual membrane with the smallest distance from a certain subproblem (in fact, different upper bound vector represents different subproblem) in N individual membranes to the subproblem. Line 8 defined the formula of calculating the distance. is matching strategy is implemented after transforming the main objective function, which means each upper bound vector is combined with the fittest subproblem throughout. By implementing this strategy, the excellent individuals will not be eliminated because of a single evaluation standard; thus, the diversity of the membrane population is maintained.

Division Operator.
In the evolutionary process, the division operator is used to generate new individual membranes based on the existing heuristic information to ensure that the membrane population always evolves in a better direction. Algorithm 7 shows the implementation process of the division operator.
One of the parent membranes used to perform division operator comes from the EM, and another comes from the whole membrane population or the neighbor membranes of the current membrane depending on probability. In the process of division operator, the decision variables shared by the two parental membranes are retained firstly, and then the decision variables unique to the two parental membranes are removed or retained by equal probability. In removing or keeping decision variables, the selection is made based on the score of the decision variable, the decision variable with a x 5 , x 7 , smaller score will be retained, and the decision variable with a larger score will be removed. Finally, to avoid the membrane population falling into the local optimum, the decision variables are released from the membrane to the environment or copied from the environment to the membrane with equal probability. Similarly, the decision variables with a larger score in the individual membrane will be released, while the decision variables with a smaller score in the environment will be copied into the individual membrane.

Fusion
Operator. When the division operator generates a new individual membrane, it may not perform well for the subproblem corresponding to the upper bound vector in its membrane, but performs well for the subproblem corresponding to another upper bound vector. So, Algorithm 8 fuses the newly generated individual membrane with the other upper bound vector, and the neighbor membrane will also be updated. Algorithm 8 contains two parts. Firstly, the newly generated individual membrane is fused with another upper bound vector suitable itself, to avoid the waste of the newly generated excellent individual membrane because it does not adapt to its own upper bound vector. Lines 1-6 define the steps to search the upper bound vector corresponding to the newly generated individual membrane, which is conducive to the convergence of the membrane population.
en, according to the objective function value of the newly generated individual membrane, the neighbor membranes are updated, the nondominated individual membrane is retained, and the dominated solution is removed.

Cytolysis and Selection Operator.
e selection and cytolysis operator updates the EM by calculating the dominant relationship between the newly generated individual membrane and the individual membrane in the EM. All the individual membranes dominated by the newly generated individual membrane in the external archive membrane will be dissolved. If no individual membrane dominates the newly generated individual membrane in the EM, the newly generated individual membrane will be added to the EM. Algorithm 9 shows the selection and cytolysis operator in detail. ese four evolutionary operators occurred during the evolution of the entire membrane population. e division operator produces better individual membrane, the fusion operator matches the individual membrane with the appropriate upper bound vector, and the selection and cytolysis operator guarantees the quality of the EM. e execution of different operators promotes the exchange of information between membranes and makes the entire membrane population continue to evolve.
In general, MEA-CNDP is a new MEA framework for solving bi-CNDP. e main features of the MEA-CNDP algorithm are as follows: (1) since each membrane represents an individual, the entire solving process can be realized entirely through membrane evolving. (2) e evolutionary process includes not only the evolution of the objects in the membrane but also the evolution of the membrane structure, so that excellent individuals can always be maintained. (3) Since the substance stored in the membrane can be any object, the complexity of encoding different evolutionary objects is avoided. (4) e evolutionary strategy guarantees the maximum theoretical parallelism of the evolutionary process and greatly improves the efficiency of the algorithm. To verify the effectiveness and efficiency of MEA-CNDP, based on four different types and sixteen different sizes of CNDP instances, this paper designs relative comparative experiments of the MEA-CNDP algorithm and some other algorithms, and the relevant results will be introduced in Section 4.

Numerical Experiments
is section is devoted to design related experiments to verify the effectiveness of the proposed MEA-CNDP algorithm and compare it with other existing algorithms on bi-CNDP. Specifically, four different types of benchmark problems and datasets of CNDP are outlined firstly. As each type of problem possesses four different sizes of datasets, 16 benchmark instances, the parameter characteristics of all instances are introduced as well. To compare the performance of different algorithms objectively, the experimental performance measures are given secondly. en, the parameters of all the algorithms involved are provided, as well as the experimental environment. Finally, the experimental results of the MEA-CNDP and other algorithms are shown, and further analysis of the results proves the effectiveness of the MEA-CNDP.

e Benchmark Problems and Dataset of CNDP.
Since there are no specific datasets for the biobjective critical node detection problem, all of the instances of this paper are composed of the graphs of several complex networks [27].
ese graphs are usually studied as single-objective critical node detection problem, and many effective results have been obtained. e datasets contain 16 undirected unweighted graphs, all of which are created by complex network generator algorithms. However, few studied deal with them as MOPs, which are also the motivations and inspirations of this paper. e datasets are divided into four types. Barabasi-Albert (BA) graphs start from a small number of nodes, enriches the network with the increase of nodes and edges over time, and finally forms a scale-free complex network graph, which is proved to be the easiest to process in related data sets. Watts-Strogatz (WS) graphs simulate a small world with a more intensive structure, that is, the diameter of the graph is small, and most nodes can access each other in relatively few hops, which is the most challenging. Edros-Renyi (ER) graphs are random, and every pair of vertices in the graph are randomly connected to form edges according to the probability, to form a completely random graph. Forest-Fire (FF) graphs simulate the model of fire spreading in a forest, different from BA, FF reproduces the heavy tailed distribution, and the network diameter decreases with time. It is worth noting that although these digitized instances cannot Computational Intelligence and Neuroscience 13 be reproduced in real networks, the real complex networks often show the combination of features in these datasets. For these critical node detection problems, it is vital to take into account the number of deleted nodes and the cost of deleting nodes at the same time. erefore, it is worthy to consider them from the perspective of MOPs. MEA has been proved to have its own unique advantages in solving combinatorial optimization problems [14][15][16]. Besides, the proposed MEA algorithm in this paper takes into account the sparsity and large-scale characteristics of the above application data sets, so it shows excellent performance.
e experimental results are detailed in Section 4.4.
In order to characterize the graph structure of the above datasets accurately, some related quantities are shown in Table 1. n is the number of nodes, m is the number of edges, 〈d〉 represents the average degree of nodes, which can be formulated as 〈d〉 � 2 · m/n, and nAP represents the number of articulation nodes. CC represents the value of the clustering coefficient, D represents the average shortest path length [28], |D 1 | represents the number of nodes having degree 1, and N(|D 1 |) represents the number of nodes whose neighboring nodes having degree 1 [36]. Among these characteristics, consider the number of articulation nodes nAP because the graph with larger nAP is easier to fragment. e clustering coefficient CC can describe the clustering tendency of nodes. e average shortest path length D represents the average distance between two randomly selected nodes in the graph.
Since these networks are not weighted in the original network, new benchmark instances are created by assigning a weight value to each node of each network, and the weight of each node is regarded as the cost of removing the node. e weight generation method comes from [31] and can be described as follows: (1) Assign weight to nodes randomly, for example, cost(v) ∈ [0.2, 3], ∀v ∈ V (2) Assign weight to nodes according to the degree, for example, cost(v) � log(d v ) + 0.5, ∀v ∈ V, where d v represents the degree of the node

IGD.
IGD can provide combined information on the convergence and diversity of the known solution. erefore, it is widely used to evaluate the approximate solution set of MOPs.
Let P * denote a subset of pareto-optimal solutions uniformly distributed along with the true PF and P be an approximation set obtained by multiobjective optimization algorithms. e IGD measures the average distance from each point in P * to its nearest point in P regarding the objective space. It is calculated as where d(x, P) represents the Euclidean distance from x to its nearest neighbor in P. Generally, a smaller IGD value indicates a better convergence and diversity of P. In MEA-CNDP, for calculating IGD, N nondominated solutions are selected from the external membrane using the crowding distance approach [41]. P * is chosen from the combination of all the solutions obtained by the experiments over 20 independent runs.

HV.
HV is used to calculate the volume covered by the obtained dominant individual in the objective space, especially when the true PF of the MOPs in the application are unknown. It is calculated as where n PF is the  [42] and random mutation [22] are adopted to generating new candidate solutions. Besides, the control parameters of the related generation operation are the same as those in [22], that is, the biased probability of crossover is set as 0.65, and the random mutation probability of each decision variable of the solution is set as 0.03.
For algorithms I-MOEA/D, I-MOEA/D-εC and MEA-CNDP, the size of the neighborhood is set as ⌊0.1 · N⌋, the probability of selecting an individual from the neighborhood is set as 0.9. For algorithms I-MOEA/D and I-MOEA/D-εC, the maximal number of replacement is set as ⌊0.01 · N⌋. e maximum number of iterations I is set as 2500, 4000, 6000, and 7500 when the number of nodes is n ≤ 500, 500 < n ≤ 1000, 1000 < n ≤ 2500, and 2500 < n ≤ 5000. For variants, only the best results among their various variants are considered in the following experiments, that is, in the mating pool selection strategy, the "-NP-EP" method is adopted, and the "-G" method is adopted in the replacement pool strategy. All of the experimental results of I-MOEA/D and I-MOEA/D-εC are from [31]. Tables 2-5 exhibit the total comparative experimental results of the three algorithms. e values in Tables 2 and 3 are the average IGD on all test instances with random weights and logarithmic weights, respectively. e values in Tables 4 and 5 are the average HV on all test instances with random weights and logarithmic weights, respectively. For each indicator, the bold data in each table are the best mean metric values of each instance. Noteworthy is that, for fair comparison, the same parameters in the three algorithms are given the same values. As for the specific parameter in MEA-CNDP algorithm, the experiments to detect the influence of parameters on the experimental results are designed and implemented, and therefore, the excellent parameters are selected. e experimental results in the tables are the results of selecting better parameters. Tables 2 and 3 show the results of the IGD metric values on all test instances with random weights and logarithmic weights. As shown in Table 2, in 16 test instances with random weights, the MEA-CNDP algorithm performs better than other algorithms in 11 instances. Further analysis can find that, in all instances of BA and FF of four different sizes, the MEA-CNDP algorithm also has significant advantages over other algorithms. Moreover, the MEA-CNDP algorithm shows poor results in the WS instance. In ER instance, the MEA-CNDP algorithm and I-MOEA/D-εC algorithm show competitive performance. e experimental results in Table 3 show that, in 16 test instances with logarithmic weights, the MEA-CNDP algorithm performs better than other algorithms in 10 instances. Specifically, in 16 instances, the MEA-CNDP algorithm shows absolute advantage in BA and FF instances with all sizes of graphs and shows slightly worse results in WS and ER instances. Tables 4 and 5 show the results of the HV metric values on all test instances with random weights and logarithmic weights. Different from the experimental results of the IGD metric, it is evident that the performance of the MEA-CNDP algorithm in the HV metric is far better than the other algorithms. Specifically, the MEA-CNDP algorithm achieves the best solution on all sizes of 16 instances with random weights and achieves the best on 15 out of 16 instances with logarithmic weights, proving that the MEA-CNDP algorithm possesses good convergence and population diversity.
It can be seen from the results of the above comparative experiments that, MEA-CNDP algorithm has good performance for most of the instances, except for the WS and ER instances.
To better measure the performance of MEA-CNDP more accurately, further experiments of I-MOEA/D-εC and MEA-CNDP algorithms are carried out. Figures 9 and 10 exhibit the degree of improvement of the MEA-CNDP algorithm relative to the I-MOEA/D-εC algorithm. In Figure 9, PIR IGD is calculated by formula (8), and IGD(A) represents the IGD metric of algorithm A. PIR IGD represents the relative change rate of the IGD metric of the algorithms MEA-CNDP and MOEA/D-εC as the number of decision variables increases. e smaller the IGD, the better the solution. So, MEA-CNDP shows worse results than I-MOEA/D-εC in smallscale instances of WS and ER. However, with the increase of the number of decision variables, the MEA-CNDP algorithm shows better results.
In Figure 10, PIR HV is calculated by formula (9), and HV(A) represents the HV metric of algorithm A. PIR HV represents the relative change rate of the HV metric. e larger the HV, the better the solution, which means that the MEA-CNDP algorithm is always better than the I-MOEA/ D-εC algorithm in HV metric, and this advantage becomes more obvious as the number of decision variables increases.
us, it can be seen that the MEA-CNDP algorithm possesses better adaptability for large-scale MOPs.
In addition, some potential knowledge can be derived from the initialization strategy of the MEA-CNDP. Algorithms 3 and 4 introduce the evaluation method of decision variables and the initial population generation strategy, respectively. Among them, for evaluating the score value of decision variables, Algorithm 3 eliminates the consideration of nodes with degree 1 and considers those with a degree greater than 1 because deleting the nodes with degree 1 in the graph will only have a small impact on calculating the pairwise connectivity of the residual graph. It can be seen from Table 1 that there are no nodes with degree 1 in WS instances and only a few nodes with degree 1 in ER instances; in other words, WS and ER instances are relatively denser graphs. Furthermore, Algorithms 3 and 4 consider the ratio of the nondominated front number to the cost of deleting nodes as the scores of decision variables, which also shows that the MEA-CNDP algorithm is more suitable for sparse graphs.
An excellent initial population will have a significant impact on the subsequent solution space search. As can be seen, the I-MOEA/D-εC algorithm uses the random initialization strategy; in contrast, MEA-CNDP designs the above specific initialization strategy. Figure 11 shows the population generated by the initialization strategy of MEA-CNDP and random initialization strategy in the objective space, and each population contains 50 solutions for the four types of instances of 500 decision variables with random weights, where f 1 and f 2 are two objective function values. Not difficult to see from Figure 11, the random population initialization strategy often has poor convergence and population diversity, so that it is easy to cause the subsequent search to fall into the local optimum, thereby deteriorating the entire algorithm process. In contrast, the population initialization strategy of the MEA-CNDP algorithm considers the characteristics of each decision variable itself. It performs the necessary screening of the solution when the population is initialized, so that the entire initial population has a higher quality, which is undoubtedly beneficial to produce a better solution. is also shows that our strategy has greater advantages when facing large-scale sparse graphs such as BA and FF instances.

Conclusion and Future Work
As a classical problem in combinatorial optimization, CNDP is often studied as a single-objective optimization problem. However, this kind of research method has some limitations and cannot cover all the problems. In addition, as a kind of new evolutionary algorithm, MEA has its unique advantages in solving such combinatorial optimization problems with graph properties. erefore, this paper is devoted to research a membrane evolution algorithm to solve the problem of bi-CNDP (MEA-CNDP). In MEA-CNDP, a population initialization method based on decision variable evaluation is proposed. en, some measures to enhance the efficiency of the MEA-CNDP algorithm are presented, including the main objective function transforming, membrane and subproblem matching strategies. Finally, according to the algorithm framework of MEA-CNDP, we design and implement the correlation division, fusion, cytolysis, and selection operators. To verify the effectiveness of MEA-CNDP, the comparative experiments are designed. e experimental results show that MEA-CNDP has good performance in solving bi-CNDP.
For future work, there are two potential research directions. At first, the population initialization strategy is worth further consideration, a more accurate method to calculate the score of the decision variable can be designed, and updating the score of the decision variable in the evolutionary process is considerable as well. Another interesting research direction is considering some local search strategies in the design of evolutionary operators, which can further enhance the quality of solutions.

Data Availability
All the original experimental data in this article are from http://individual.utoronto.ca/mventresca/cnd.html.

Conflicts of Interest
e authors declare that they have no conflicts of interest.