Weights Adaptation Optimization of Heterogeneous Epidemic Spreading Networks: A Constrained Cooperative Coevolution Strategy

1 W eights Adaptation Optimization of Heterogeneous Epidemic Spreading Networks: A Constrained Cooperati v e Coe volution Strate gy Y un Feng and Bing-Chuan W ang Abstract —In this paper , the dynamic constrained optimiza- tion pr oblem of weights adaptation for heter ogeneous epidemic spreading networks is in vestigated. Due to the powerful ability of searching global optimum, e volutionary algorithms are employed as the optimizers. One major difficulty is that the dimension of the problem is incr easing exponentially with the network size and most existing ev olutionary algorithms cannot achie ve satisfiable performance on large-scale optimization problems. T o address this issue, a novel constrained cooperative coevolution ( C 3 ) strategy , which can separate the original lar ge-scale problem into differ ent subcomponents, is employed to achieve the trade- off between the constraint and objective function. Index T erms —Evolutionary Computation, Constrained Opti- mization, Epidemic Spreading, W eights Adaptation. I . I N T R O D U C T I O N E PIDEMIC spreading ov er complex netw orks [1] has attracted lots of attention since the pioneering work [2] of Daniel Bernoulli in 1760. Many researches hav e been focused on the mathematical modeling of disease spreading process, classical epidemic models such as the susceptible-infected- susceptible (SIS) [3], [4], [5], [6], and the susceptible-infected- recov ered (SIR) [7] model have been well studied for decades. Since the spreading of disease (such as AIDS, SARS, etc.) may cause numerous damages to the human society , dev eloping control policies for epidemic spreading process is of great significance, with potential applications in public health. As pointed out in [8], the two common strategies to suppress epidemic spreading scale are increasing the recovery rate and decreasing the infection rate. For example, in [9], the PID control laws were implemented for the classical SIR model where the v accination rate is the control v ariable. Despite these innov ative results, these studies did not take the “budget” or the so-called “control cost” into consideration, which must be taken into account in real-world scenarios . In recent years, optimal control of epidemic spreading [10], [11], [12], [13] has gained more and more attentions. For example, in [14], an optimal control strategy was designed for the vaccina- tion, quarantine and treatment actions. F or inhomogeneous epidemic dynamics, the optimal control problem was studied in [15]. Y un Feng is with the Department of Systems Engineering and Engineering Management, City Uni versity of Hong Kong, Ko wloon, Hong Kong, and also with the Central South Uni versity , Changsha, 410083, P . R. China. E-mail: yun.feng@my .cityu.edu.hk Bing-Chuan W ang is with the School of Automation, Central South Univer - sity , Changsha, 410083, P . R. China. E-mail: wangbingchuancityu@gmail.com (Corresponding author: BC W ang) Besides the above control strategies which are focusing on the epidemic spreading parameters, recently another control strategy which aims at adjusting topology of the underlying network was in vestigated [16]. Acting as the “bridge” for epidemic spreading from infected individuals to healthy ones, the network topology determines the epidemic transmission efficienc y . A quantitati ve parameter that defines the strength of two nodes in a network is the v alue of weight between them. Since the physical meaning of the weight can be described by the contact frequency of two indi viduals, the intuitiv e idea of controlling the weights is more natural and practical than controlling the spreading parameters. In [17], an individual-based weight adaptation mechanism in which individuals’ contact strength is adaptable depending on the lev el of contagion spreading over the network was proposed. An optimal control formulation was also presented to address the trade-off between the global infected lev el and the local weight adaptation cost corresponding to the topology of the underlying contact network. In [16] and [17], the objective function contains both the infection cost and control cost. Howe ver , the problem of minimizing the infection cost with giv en fix ed number of budget, i.e., the control cost, has not been studied, which is more practical than the unconstrained optimization problems. This moti v ates us to in vestigate the constrained optimization problem of weights adaptation. Meanwhile, when solving the optimal control problems in both [16] and [17], the forward-backward sweep method (FBSM) [18] was used to find the numerical solutions, which is an indirect method itself. The intuitiv e idea of FBSM is that the initial value problem of the state equations is solved forward in time, using an estimate for the control and adjoint variables. Then the adjoint final value problem is solved backwards in time. These complicated procedures add sev eral difficulties to the optimal problem, which can be summarized as follows: • Need to compute various partial deriv ati ves of the Hamil- tonian and solve additional differential equations, which introduce more errors for the optimization problem. • Need to make an initial guess of the adjoint v ariables, the sensitivity of the method to changes in initial guesses. In [18], the authors found that the FBSM method fails to terminate under some circumstances. Evolutionary algorithms (EAs) [19], [20] which are inspired by nature, ha ve sho wn po werful searching ability for the global optimum with fe w restrictions. Also, since EAs are direct 2 methods, the implementations are much simpler than other con ventional methods. EAs have been widely used in the community of network science and engineering [21], [22], [23]. For example, a nov el memetic algorithm which can preserve the community structure is proposed for the network embedding problem in [21]. As one of the most po werful ev olutionary algorithm, differential e volution (DE) [24], [25], [26] has sho wn superior performance ov er other heuristic algo- rithms on very complex searching and optimization problems. The control parameters of DE are few and it is highly efficient [24]. This motiv ates us to solve the constrained optimization problem based on DE. Howe ver , since the weights adaptation in volves all links in the network for a given period of time [17], the dimension of the optimal solution is relativ ely high. Solving this kind of large-scale optimization problem [27] is a challenging problem in the community of ev olutionary computation. In addition, how to achie ve the trade-off between the constraint and objective function is another difficulty . Motiv ated by the above considerations, a dynamic con- strained optimization problem of weights adaptation for het- erogeneous epidemic spreading networks based on SIS model (In this manuscript, we only consider SIS model and SIR model will be studied in the future work) is formulated. T o deal with the high-dimensional optimization problem, a nov el constrained cooperativ e coev olution ( C 3 ) strategy which can separate the original high-dimensional search space into some lo w-dimensional ones by random grouping strategy is proposed. The  constraint-handling technique is employed to achiev e the trade-off between the constraint and objecti ve function. Moreover , as a commonly used variant of the clas- sical DE, the differential ev olution with neighborhood search (NSDE) [28] algorithm is employed as the optimizer for these sub-problems. The main contrib utions of this paper can be summarized as follows: • A dynamic constrained optimization problem of weights adaptation for heterogeneous epidemic spreading net- works is formulated. • Evolutionary computation techniques with strong search- ing ability as well as quite a small amount of demands are applied to solve this kind of problem. • A nov el constrained cooperativ e coev olution ( C 3 ) strat- egy is tailored for this real-world large-scale optimization problem. • It is easy to implement and the optimization process can be done once the basic network and epidemic param- eters are given. For the feedback control strategy [29], information about the number of susceptible or infected individuals is needed for the control law updating while the proposed method does not require such information. The proposed strategy has strong potentials to be applied to real-world scenarios for disease control. Actually the most widely and successfully used strategy in real-world, quaran- tine, is a special form of weights adaptation. In addition, the proposed strategy has no restrictions on the spreading network. The rest of this paper is organized as follows. The model description and problem formulation are given in Section II. Differential ev olution is briefly introduced in Section III. In Section IV, the methodology is giv en in detail. Some numerical experiments are presented in Section V. Finally , this paper is concluded in Section VI. I I . M O D E L D E S C R I P T I O N A N D P RO B L E M F O R M U L A T I O N A. Heter ogeneous W eighted SIS-based Network Model In the heterogeneous weighted SIS-based network, all nodes can be classified into two possible states according to their health status, that is, susceptible and infected. Susceptible individuals can be infected by infected individuals through the links between them and in state X i ( t ) = 1 . Meanwhile, infected indi viduals can be cured and become susceptible ag ain and in state X i ( t ) = 0 . T o be more specific, ev ery node i at time t is infected with probability Pr [ X i ( t ) = 1] and susceptible with probability Pr [ X i ( t ) = 0] . At each time t , a node can only be in either of these two states, thus Pr [ X i ( t ) = 1]+ Pr [ X i ( t ) = 0] = 1 . The state for each individual at time t is independent. Then the following heterogeneous weighted SIS-based net- work model is obtained from the N-intertwined mean-field approximation (NIMF A) [30], [31], [17]: ˙ p i ( t ) = − γ i p i ( t )+(1 − p i ( t )) N X j =1 w ij ( t ) β j p j ( t ) , i = 1 , · · · , N (1) where p i ( t ) ∈ [0 , 1] denotes the probability of node i being infected at time t ≥ 0 . The infection and curing rates β i ≥ 0 and γ i ≥ 0 for each node i in the network are described by two independent Poisson processes; w ij ( t ) ∈ [0 , 1] denotes the weight of edge from node j to node i ; N is the number of nodes in the network. Then (1) can be rewritten in the following compact form: ˙ p ( t ) = ( W ( t ) B − D ) p ( t ) − P ( t ) W ( t ) B p ( t ) , (2) where p ( t ) = [ p 1 ( t ) , · · · , p N ( t )] T , W ( t ) = [ w ij ( t )] N × N , B = diag [ β 1 , · · · , β N ] , D = diag [ γ 1 , · · · , γ N ] and P ( t ) = diag [ p 1 ( t ) , · · · , p N ( t )] . Here “ diag ” denotes a diagonal matrix. T o be noticed, the network considered here is a directed one and the weights of which are in association with the infected lev el in the network as it is shown in (2). Define w ii = 0 for all node i so that self-loop is not considered here. B. Pr oblem F ormulation As it is shown in (2), the infected lev el can be controlled by the weights adaptation in the network. Howe ver , bearing the cost of weights adaptation in mind, a natural dynamic constrained optimization problem is formed as follows: min w ∈ W f = Z T 0 { N X i =1 f i ( p i ( t )) } dt (3) 3 s.t. : g = Z T 0 { N X i =1 N X j =1 g ij ( w ij ( t ) − w 0 ij ) } dt − C ≤ 0; ˙ p i ( t ) = (1 − p i ( t )) N X j =1 w ij ( t ) β j p j ( t ) − γ i p i ( t ); p i (0) = p 0 ( i ) , 0 ≤ p i ( t ) ≤ 1 , 0 ≤ w ij ( t ) ≤ 1 , 1 ≤ i ≤ N , 1 ≤ j ≤ N . where W is the set of all admissible weights; f i ( p i ( t )) denotes the infection cost function for each indi vidual i , the objectiv e function f denotes the total infection cost for all individuals in the considered period [0 , T ] ; g ij ( w ij ( t ) − w 0 ij ) denotes the cost function for weights adaptation and w 0 ij denotes the initial weight at t = 0 , C is a constant that characterizes the maximum cost of weight adaptation, g is the inequality constraint. The dynamic constrained optimization problem is inter- preted as: how to design the adaptiv e weights { w ij ( t ) } for all links from t = 0 to t = T such that the infected level in the network can be suppressed to the maximally extent under giv en budgets. I I I . D I FF E R E N T I A L E V O L U T IO N Differential ev olution (DE) [32], [33] is used as a base optimizer to solve the dynamic constrained optimization problem in (3), which is arguably one of the most powerful stochastic real-parameter optimization algorithms in current use. Different from traditional EAs, DE uses difference of individual trial solutions to explore the objective function landscape. T ypically , DE consists of four stages: initialization, mutation, crossover , and selection. Step 1: Initialization The ultimate goal of DE is to find a global optimum point in a D -dimensional real parameter space R D . In the beginning, a randomly initiated population P which consists of N P (size of population) D -dimensional individuals is generated. N P denotes the size of population where each individual with size N × ( N − 1) × ( T − 1) denotes a solution that consists of the weighs of ev ery directed link in the spreading network from time t = 1 to time T − 1 . The subsequent generations in DE are denoted by G = 1 , · · · , G max . The population at generation G is denoted as follows: P G = { ~ x 1 ,G , · · · , ~ x N P ,G } , where ~ x i,G is the i th target vector in current generation: ~ x i,G = [ x 1 ,i,G , · · · , x D,i,G ] . Initially ( G = 0 ), the population should be uniformly random- ized within the search space constrained by the maximum and minimum bounds: ~ x min = [ x 1 , min , · · · , x D, min ] and ~ x max = [ x 1 , max , · · · , x D, max ] . Hence the initialization is formulated as follows: x j,i, 0 = x j, min + r and i,j [0 , 1] · ( x j, max − x j, min ) , where r and i,j [0 , 1] is a uniformly distributed number between 0 and 1 . Step 2: Mutation The mutation operator aims to create a mutant vector for each target vector through utilizing the differential information of pairwise individuals. The following mutation operator is adopted in this paper . DE/current-to-best/1: ~ v i,G = ~ x i,G + F ( ~ x best,G − ~ x i,G ) + F ( ~ x r i 1 ,G − ~ x r i 2 ,G ) , where i = 1 , · · · , N P , ~ v i,G = [ v 1 ,i,G , · · · , v D,i,G ] is the mutant vector , r i 1 and r i 2 are mutually exclusi ve integers randomly chosen from [1 , N P ] \ i , ~ x best,G is the best individual in the current population, and the scaling factor F is a positiv e control parameter for scaling the difference vectors. Step 3: Cr ossover The crossover operator is employed to enhance the diversity of the population. The trial vector ~ u i,G = [ u 1 ,i,G , · · · , u D,i,G ] is formed by exchanging components between the target vector ~ x i,G and the mutant vector ~ v i,G . The follo wing binomial crossov er is utilized: u j,i,G =  v j,i,G , if ( r and i,j [0 , 1] ≤ C r or j = j rand ) x j,i,G , otherwise. where i = 1 , · · · , N P , j = 1 , · · · , D , j rand is a randomly chosen index from [1 , D ] . C r is the crossover rate. Step 4: Selection In the selection step, the target vector ~ x i,G is compared with the trial vector ~ u i,G , the better one can be survived in the next generation ~ x i,G +1 =  ~ u i,G , if f ( ~ u i,G ) < f ( ~ x i,G ) ~ x i,G , otherwise. I V . M E T H O D O L O G Y A. Encoding Mechanism As described above, the dimension of the decision space D is related to the number of links in the network and the time length T . T o be more specific, for a directed network with N nodes, we hav e D = N × ( N − 1) × ( T − 1) . (4) T o better illustrate the encoding mechanism, a schematic graph is presented in Figure 1. The network consists of N = 4 nodes, where blue and red nodes represent susceptible and infected individuals, respecti vely . ~ x i,G ( t ) denotes the i th target vector in the G -th generation, with a dimension of N × ( N − 1) . Each element of ~ x i,G ( t ) corresponds to a specific element of the weight matrix W i,G ( t ) at time t in the manner presented in Figure 1. Bearing the infection cost and the weights adaption constraint in mind, the weight matrix W i,G ( t + 1) at time t + 1 is ev olved adaptiv ely to balance the objectiv e function 4 , 0.1 0.3 0.2 0.6 0.5 0.5 0. 0 0 0 ( ) 7 0.3 0.4 0.8 0. .1 0 9 0 i G t W æ ö ç ÷ ç ÷ = ç ÷ ç ÷ è ø 1 2 3 4 1 2 3 4 , 0.3 0.2 0.5 0.4 0.4 0.3 0 0 0 0 .3 0.2 0.5 0.6 0.7 0.3 ( ) 0 1 i G W t æ ö ç ÷ ç ÷ = ç ç è ø + ÷ ÷ , 0.1 0.3 0.2 0.6 0.5 0.5 0.7 0.3 0.4 0.8 0.9 0. ( [ ) ] 1 i G x t = r , 0.3 0.2 0.5 0.4 ( 1 ) 0.4 0.3 0.3 0.2 0.5 0.6 0.7 [ ] 0.3 i G x t = + r , , , , , ( ) ( 1 ) ( 1 ) ( 1 ) i G i G G i G i i G x x T x t x t x æ ö ç ÷ ç ÷ ç ÷ ç ÷ = ç ÷ ç ÷ ç ÷ ç ÷ - è ø + r M M r r r r , , 1, G i G NP G G x P x x æ ö ç ÷ ç ÷ ç ÷ = ç ÷ ç ÷ ç ÷ è ø r M M r r Fig. 1. Schematic of the encoding mechanism. and constraint. Considering the time length T , the i th target vector at G -th generation is formulated as: ~ x i,G = [ ~ x i,G (1) , · · · , ~ x i,G ( t ) , ~ x i,G ( t + 1) , · · · , ~ x i,G ( T − 1)] . T o be noticed, the health status of all individuals in the net- work is varying from time to time. For example, as it is shown in Figure 1, node 1 is infected at time t while it is cured and become susceptible at time t + 1 . Therefore, this is a dynamic optimization problem considering the interactions between the epidemic spreading process and weights adaptation. Remark IV .1. W ith this encoding mechanism, one inevitable pr oblem is that the dimension of the decision space increases exponentially with the size of the network N as it is shown in (4). Ther efore, this constrained optimization pr oblem (3) may suffer s from the “curse of dimensionality”, which implies that most of the EAs’ performance deteriorates rapidly as the incr easing of the dimensionality of the searc h space. B. Differ ential Evolution with Neighborhood Sear ch (NSDE) As a v ariant of classical DE introduced in Section III, NSDE [28] is effecti ve in escaping from local optima when searching in circumstances without knowing the preferred step size. The main difference between NSDE and DE is that the neighborhood search (NS) strategy is utilized, which is a typical technique in ev olutionary programming (EP) [34]. T o be more specific, the scaling factor F in classical DE is replaced in the following manner: F i =  N i (0 . 5 , 0 . 5) , if ( r and i [0 , 1] < f p ) δ i , otherwise. where N i (0 . 5 , 0 . 5) is a Gaussian random number with mean 0.5 and standard de viation 0.5, and δ i is a Cauchy random variable with scale parameter t = 1 . In NSDE, the parameter f p was set to a constant number 0.5. C.  Constraint-handling T echnique The  constraint-handling technique [35], which is adopted by the winner of IEEE CEC2010 competition, is utilized to compare two target vectors ~ x i and ~ x j . T o be more specific, ~ x i is better than ~ x j if the following conditions are satisfied:      f ( ~ x i ) < f ( ~ x j ) , if G ( ~ x i ) ≤  ∧ G ( ~ x j ) ≤  (5a) f ( ~ x i ) < f ( ~ x j ) , if G ( ~ x i ) = G ( ~ x j ) (5b) G ( ~ x i ) < G ( ~ x j ) , otherwise. (5c) where G ( ~ x ) denotes the degree of constraint violation on the constraint as follows: G ( ~ x ) = max(0 , g ( ~ x )) . In Eq. 5a,  is designed to decrease with the increasing of the generation G as follows [36]:  =     0 (1 − G G max ) cp , if G ≤ Gc 0 , otherwise. cp = − log  0 + λ log(1 − Gc G max ) , where  0 is the maximal degree of constraint violation of the initial population; Gc is a parameter to truncate the value of  ; λ is set to be 10 in this paper . D. A Constrained Cooperative Coevolution ( C 3 ) Strate gy In this subsection, a novel coev olution strategy is dev el- oped to solve this dynamic constrained optimization problem motiv ated by the DECC-G algorithm in [27]. The core of this coev olution strate gy is “divide-and-conquer”. That is, the original D -dimensional search space is divided into N s number of D s -dimensional ones by random grouping strategy , where N s × D s = D . 5 The constrained cooperative coevolution ( C 3 ) framew ork for this problem can be summarized as follows: (1) Set i = 1 to start a new cycle . (2) Decompose an original D -dimensional target vector into N s low-dimensional subcomponents with dimension D s randomly , i.e. D = N s × D s . Here “randomly” indicates that each dimension in the original target vector has the same probability to be assigned into any of the N s subcomponents. (3) Optimize the j th subcomponent with NSDE and  constraint-handling technique introduced in Section IV -B and Section IV -C for a predefined number of fitness e valuations (FEs). (4) if j < N s then j + + and go to Step 3. (5) Stop if halting criteria are satisfied; otherwise go to Step (1) for the next cycle . Remark IV .2. The pr obability of C 3 strate gy to assign two interacting variables x i and x j into a single subcomponent for at least k cycles is: P k = K X l = k  K l  1 N s  l  1 − 1 N s  K − l wher e K is the total number of cycles and N s is the number of subcomponents. In each separate cycle, the probability to assign two inter- acting variables x i and x j into a single subcomponent is p =  N s 1  N 2 s = 1 N s . Let p k denotes the pr obability to assign x i and x j into a single subcomponent for exactly k cycles. Obviously , p k satisfies the binomial distribution, so: p k =  K l  p l (1 − p ) K − l =  K l  1 N s  l  1 − 1 N s  K − l . Thus, P k = K X l = k p k = K X l = k  K l  1 N s  l  1 − 1 N s  K − l . Giv en network size N = 20 and T = 10 , then D = N × ( N − 1) × ( T − 1) = 3420 . Select D s = N × ( N − 1) = 380 , then the number of subcomponents N s = 9 . When the number of cycles K = 50 , we have: P 1 = 1 − (1 − 1 9 ) 50 = 0 . 9972 . P 2 = P 1 −  50 1  × 1 9 × (1 − 1 9 ) 49 = 0 . 9799 . These results demonstrate that the C 3 strategy has relatively high probabilities to optimize interacting v ariables in a single subcomponent for at least one or two cycles. E. NSDE under the C 3 frame work Considering NSDE as the base optimizer for the subcom- ponents, it is straightforward to obtain the NSDE with con- strained cooperative coev olution algorithm, denoted as NSDE- C 3 . The pesudocode of NSDE- C 3 is given in Algorithm 1. Algorithm 1: NSDE- C 3 1 /* Initialization */ 2 pop (1 : N P , 1 : D ) ← r and ( popsiz e, D ) ; 3 ( best, bestv al ) ← ev aluate ( pop ) ; 4 for i = 1 : cy cl es do 5 index (1 : D ) ← r andper m ( D ) ; 6 for j = 1 : N s do 7 k ← ( j − 1) × D s + 1 ; 8 l ← j × D s ; 9 subpop (1 : N P, 1 : D s ) ← pop (: , index ( k : l )) ; 10 /* Use sub-optimizer */ 11 subpop ← N S DE ( best, subpop, F E s ) ; 12 pop (: , index ( k : l )) ← subpop (1 : N P , 1 : D s ) ; 13 ( best, bestv al ) ← ev aluate ( pop ) ; 14 return pop ( best, :) , bestv al V . N U M E R I C A L E X P E R I M E N T S T o validate the effectiv eness of the proposed method, a synthetic network is obtained as the underlying netw ork. When solving the dynamic constrained optimization problem in (3), each optimization algorithm is used for 25 independent runs. A. Datasets For validation, one of the most famous synthetic complex network-the B ´ arabasi-Albert (BA) network [37] is used. When constructing the BA network, initially m 0 = 5 fully connected nodes are placed in the network. A new node is connected to m = m 0 = 5 existing nodes with probability proportional to the degree of them at each step [38]. A total number of N = 20 nodes is fixed for the B A network in the simulations, originally this B A network is a bidirectional network. The topologies of the B A network is presented in Figure 2, where circles represent nodes in the network and their sizes are proportional to degrees. T able I belo w summarizes the topological features of the network, where N , < k >, < C >, d denote the total number of nodes, av erage degree, av erage clustering coef ficient and density of the network, respectiv ely . T ABLE I T OP O L OG I C A L F E A T U R ES O F T H E A N ALY ZE D N E TW O R K . Network name N < k > < C > d B ´ arabasi-Albert 20 8.5 0.519 0.447 B. P arameters Selection For all experiments, the objective function for individual i at time t adopts the form: f i ( p i ( t )) = p p i ( t ) . Meanwhile, the cost function for weights adaptation for a pair of individuals i and j is: g ij ( w ij ( t )) = ( w ij ( t ) − w 0 ij ) 2 . 6 The value of the constraint C is selected as 700. The objective function f i ( p i ( t )) and the cost function g ij ( w ij ( t )) adopt the current form for simple illustration. They represent one possible situation, there exists other sets of f i ( p i ( t )) and g ij ( w ij ( t )) in real-world scenarios. Since the weights adaptation optimization problem is still an open- problem and there lacks standard definition of both the ob- jectiv e function and the cost function, we adopted the current form similar to the way that was done in [16] and [17]. As for the value of the constraint C , it was selected after the spreading network and epidemic parameters have been set. Based on the underlying network introduced in the previous subsection, some numerical simulations are conducted to illus- trate the effecti veness of the proposed method (NSDE- C 3 ). The epidemic parameters used in the simulations are pre- sented in T able II, where p i (0) , β i , γ i , T denote the initial infection state, infection rate, curing rate, and terminal time, respectiv ely . These parameters are carefully chosen to make sure that the infected le vel is relatively high when there is no weights adaptation, hence the effects of weights adaptation on epidemic spreading can be easily observed. Remark V .1. According to Theor em 2 in [39], assume that p (0) 6 = 0 , then the metastable state p ? is globally asymptoti- cally stable when λ max ( AB − D ) > 0 wher e λ max ( AB − D ) is the larg est eigen value-the spectral radius of the matrix AB − D . Referring to the epidemic parameters in T able II, it is easily obtained that λ max ( AB − D ) = 3 . 5167 > 0 for the constructed BA network, hence the infected level is r elatively high without weights adaptation. Regarding to the algorithm parameters, the population size N P , the number of FEs and the crossover rate for both NSDE and NSDE- C 3 algorithms are the same. For NSDE- C 3 strategy , the dimension of the subcomponents D s is selected as N ( N − 1) . Fig. 2. B ´ arabasi-Albert network T ABLE II E P ID E M I C A N D A L GO R I TH M PA RA M E TE R S U S E D I N T H E S I M U LAT IO N S . Network name p i (0) β i γ i T N P FEs C r B ´ arabasi-Albert 0.153 0.4 0.3 10 350 6.30e+06 0.9 C. Results For B ´ arabasi-Albert network, the simulation results are presented in Figure 3(a), 3(b), 3(c) and 3(d). Both the pro- posed NSDE- C 3 and the existing NSDE method are used. Meanwhile, two control groups are added for comparison, i.e. “No adaptation” and “constant adaptation” strategy . The formulation of these two strategies at time t are as follows: No adaptation: w ij ( t ) = w 0 ij ; Constant adaptation: w ij ( t ) = c × w 0 ij , where g = Z T 0 { N X i =1 N X j =1 g ij ( w ij ( t ) − w 0 ij ) } dt − C = 0 . In this manner , these two strategies are employed as the baselines for comparison. “No adaptation” means that the weights remain unchanged in the considered time period, while “Constant adaptation” refers to a “discount” on the original weights w 0 ij with the budget C fully used. For the epidemic parameters in T able II and C = 700 , the constant adaptation ratio can be calculated as c = 0 . 33 . Figure 3(a) and 3(b) illustrate the ev olution process of the mean best fitness v alue and constraint violation with respect to generations over 25 independent runs, respectiv ely . For the “No adaptation” and “Constant adaptation” strategy , the solutions are determined initially . Hence the value of the objectiv e function and constraint violation for these two strategies remain unchanged. For the “No adaptation” and “Constant adaptation” strategies, the constraint violation re- mains zero. Combing Figure 3(a) and 3(b), it is evident that the proposed NSDE- C 3 outperforms NSDE on B A network. Moreov er , NSDE- C 3 outperforms both the “No adaptation” and “Constant adaptation” strategies. Remark V .2. One inter esting phenomenon to be noticed in F igure 3(b) is that the evolution of constraint violation for NSDE- C 3 almost coincides with that of NSDE. However , the e volution of the mean best fitness value for these two algorithms ar e quite differ ent. This is due to the use of  constraint-handling technique in Section IV -C, which puts mor e emphasis on the constraint violation than the fitness value for earlier generations. Figure 3(c) and 3(d) demonstrate the ev olution process of two main indices, i.e. infected le vel and total weights over time, which characterize the epidemic spreading and weights adaptation lev el, respecti vely . The definition of theses two 7 0 3 0 0 0 6 0 0 0 9 0 0 0 1 2 0 0 0 1 5 0 0 0 1 8 0 0 0 1 2 0 1 5 0 M e a n b e s t f i t n e s s v a l u e G e n e r a t i o n s N o a d a p t a t i o n N S D E C o n s t a n t a d a p t a t i o n N S D E - C 3 (a) Fitness value 0 3 0 0 0 6 0 0 0 9 0 0 0 1 2 0 0 0 1 5 0 0 0 1 8 0 0 0 0 2 0 0 4 0 0 0 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 2 0 0 4 0 0 N o a d a p t a t i o n N S D E C o n s t a n t a d a p t a t i o n N S D E - C 3 M e a n b e s t c o n s t r a i n t v i o l a t i o n G e n e r a t i o n s (b) Constraint violation 0 2 4 6 8 1 0 0 . 2 0 . 4 0 . 6 0 . 8 1 . 0 N o a d a p t a t i o n N S D E C o n s t a n t a d a p t a t i o n N S D E - C 3 I n f e c t e d L e v e l T i m e (c) Epidemic ev olution 0 2 4 6 8 1 0 4 0 8 0 1 2 0 1 6 0 2 0 0 N o a d a p t a t i o n N S D E C o n s t a n t a d a p t a t i o n N S D E - C 3 T o t a l w e i g h t s T i m e (d) W eights adapation Fig. 3. Simulation results on B A network for 25 independent runs indices are as follows: I ( t ) = E ( p i ( t )) = 1 N N X i =1 p i ( t ); W ( t ) = N X i =1 N X j =1 ,j 6 = i w ij ( t ) . The physical meaning of the infected level I ( t ) and the total weights W ( t ) are the expectation of infection probability and total number of weights, respectively . In this manner , I ( t ) can be used as an index which reflects the current infected lev el among the whole population of individuals in the network. Similarly , W ( t ) reflects the weights variation of the network concerned. For infected level, “No adaptation” strategy achiev es a relatively high level, which is used as a baseline for comparison as it is demonstrated in Section V -B and Remark V .1. Ho we ver , the “Constant adaptation” strategy achiev es a lo wer infected lev el than the “No adaptation” strategy , which is mainly due to that the “Constant adaptation” strategy makes full use of the budget as it is shown in Figure 3(d). Referring to the total weights W ( t ) variation in Figure 3(d), it is obvious that W ( t ) remains unchanged at ev ery time t for the “No adaptation” strategy . For the “Constant adaptation” strategy , the value of the total weights W decreases to a lower value and keeps unchanged for the whole time period. Howe ver , for both the NSDE and NSDE- C 3 strategies, the value of total weights W ( t ) experience a decreasing process first and restore to the initial value of time t = 1 , which coincide with the results in [17]. The restoring process of W ( t ) for NSDE- C 3 is more moderate than that of NSDE, hence the infected level I ( t ) is lower at most time. One interesting observation is that the v alue of the total weights for NSDE- C 3 drops to a lower value than that of the “Constant adaptation” strategy , and restores to the initial value at time T , the infected lev els indicate that this kind of “dynamically adaptation” is more effecti ve in controlling the epidemic spreading scale. In addition, to test statistical significance, the multi-problem W ilcoxon’ s test [36] are implemented to compare these meth- ods, the results are shown in T able III. Once more, the results indicate that the proposed NSDE- C 3 is superior to the other competitors. T ABLE III W I LC OX O N ’ S T E ST O N B A N E T WO R K F O R 2 5 I N D EP E N DE N T R UN S . Algorithm Mean OFV ± Std Dev p -value NSDE- C 3 106.4530 ± 0.1506 - NSDE 127.1566 ± 1.3425 1.4157e-09 No adaptation 160.5280 ± 0.00 9.7285e-11 Constant adaptation 126.8617 ± 0.00 9.7285e-11 Remark V .3. T o be noticed, the SIS model has been analyti- cally and exactly solved in [40], [41], wher e the infection and r ecovery rate are time-varying. However , the NIMF A model is employed as the epidemic spr eading model for its simplicity . In the futur e, we may extend our results with the help of the Lie algebra method. V I . C O N C L U S I O N In this paper , a dynamic constrained optimization problem of weights adaptation for heterogeneous epidemic spreading networks is formulated. Combining constrained cooperative coev olution ( C 3 ) strategy with NSDE, a novel NSDE- C 3 algorithm is tailored for this problem. Numerical experiments 8 on a BA network showed the effecti veness of this algorithm. In future research, how to expand this algorithm for on-line implementation is an interesting yet challenging topic. More- ov er , the extension to the Mixed SI(R) epidemic dynamics in random graphs with general degree distributions and also the feedback control optimization are also our interests. Some insights on the capability of our approach are pro- vided as follows: • First, in our problem setting, the dimension of a solution to be optimized increases exponentially with the increas- ing of the network size N , which is ob vious. And the curse of dimensionality problem is inevitable under such situations. Ho wever , e ven in the e volutionary computation community , state-of-the-art methods to deal with such large scale optimization problem is questionable. (The dimension of a solution in [27], [42] is only up to 1000 while it is 3420 in the numerical simulations studied in this manuscript.) • The goal of this manuscript is not just to purse as high dimension as possible. More importantly , we focus on introducing ev olutionary algorithms into the network optimization problems and shed some light on solving such large-scale optimization problems using the idea of cooperativ e coev olution. In the numerical simulations, the number N = 20 is selected due to the limitations on the hardware. Moreover , we do consider the capacity of the proposed algorithm seriously and there are still ways to improv e. For example, for very large N , we may consider di viding the solution of each subcomponent with dimension N × ( N − 1) into small sub-subcomponents again and then deal with it in the same manner . R E F E R E N C E S [1] R. Pastor-Satorras, C. Castellano, P . V an Mieghem, and A. V espignani, “Epidemic processes in complex networks, ” Reviews of modern physics , vol. 87, no. 3, p. 925, 2015. [2] D. Bernoulli, “Essai d’une nouvelle analyse de la mortalit ´ e caus ´ ee par la petite v ´ erole, et des av antages de l’inoculation pour la pr ´ evenir , ” Histoire de l’Acad., Roy . Sci.(P aris) avec Mem , pp. 1–45, 1760. [3] Y . Feng, Q. Fan, L. Ma, and L. Ding, “Epidemic spreading on uniform networks with two interacting diseases, ” Physica A: Statistical Mechan- ics and its Applications , vol. 393, pp. 277–285, 2014. [4] B. Qu and H. W ang, “Sis epidemic spreading with heterogeneous in- fection rates, ” IEEE T ransactions on Network Science and Engineering , vol. 4, no. 3, pp. 177–186, 2017. [5] Y . Feng, L. Ding, Y .-H. Huang, and Z.-H. Guan, “Epidemic spreading on random surfer networks with infected avoidance strategy , ” Chinese Physics B , vol. 25, no. 12, p. 128903, 2016. [6] Y . Huang, L. Ding, Y . Feng, and J. Pan, “Epidemic spreading in random walkers with heterogeneous interaction radius, ” Journal of Statistical Mechanics: Theory and Experiment , vol. 2016, no. 10, p. 103501, 2016. [7] M. Nadini, A. Rizzo, and M. Porfiri, “Epidemic spreading in temporal and adaptive networks with static backbone, ” IEEE T ransactions on Network Science and Engineering , 2018. [8] C. Nowzari, V . M. Preciado, and G. J. Pappas, “ Analysis and control of epidemics: A survey of spreading processes on complex networks, ” IEEE Contr ol Systems , vol. 36, no. 1, pp. 26–46, 2016. [9] L. L. Ghezzi and C. Piccardi, “Pid control of a chaotic system: An application to an epidemiological model, ” Automatica , vol. 33, no. 2, pp. 181–191, 1997. [10] X.-J. Li, C. Li, and X. Li, “Minimizing social cost of v accinating network sis epidemics, ” IEEE T ransactions on Network Science and Engineering , no. 1, pp. 1–1, 2017. [11] C. Nowzari, V . M. Preciado, and G. J. Pappas, “Optimal resource allo- cation for control of networked epidemic models, ” IEEE T ransactions on Contr ol of Network Systems , vol. 4, no. 2, pp. 159–169, 2017. [12] H. Chen, G. Li, H. Zhang, and Z. Hou, “Optimal allocation of resources for suppressing epidemic spreading on networks, ” Physical Review E , vol. 96, no. 1, p. 012321, 2017. [13] C. Pizzuti and A. Socievole, “ A genetic algorithm for finding an optimal curing strategy for epidemic spreading in weighted networks, ” in Pr oceedings of the Genetic and Evolutionary Computation Conference . A CM, 2018, pp. 498–504. [14] D. Xu, X. Xu, Y . Xie, and C. Y ang, “Optimal control of an si vrs epidemic spreading model with virus variation based on complex networks, ” Com- munications in Nonlinear Science and Numerical Simulation , vol. 48, pp. 200–210, 2017. [15] Y . Shang, “Optimal control strategies for virus spreading in inhomoge- neous epidemic dynamics, ” Canadian Mathematical Bulletin , vol. 56, no. 3, pp. 621–629, 2013. [16] Y . Feng, L. Ding, and P . Hu, “Epidemic spreading on random surfer networks with optimal interaction radius, ” Communications in Nonlinear Science and Numerical Simulation , vol. 56, pp. 344–353, 2018. [17] P . Hu, L. Ding, and T . Hadzibegano vic, “Indi vidual-based optimal weight adaptation for heterogeneous epidemic spreading networks, ” Communications in Nonlinear Science and Numerical Simulation , vol. 63, pp. 339–355, 2018. [18] M. McAsey , L. Mou, and W . Han, “Con ver gence of the forward- backward sweep method in optimal control, ” Computational Optimiza- tion and Applications , vol. 53, no. 1, pp. 207–226, 2012. [19] T . B ¨ ack, D. B. Fogel, and Z. Michalewicz, Handbook of evolutionary computation . CRC Press, 1997. [20] M. ˇ Crepin ˇ sek, S.-H. Liu, and M. Mernik, “Exploration and exploitation in e volutionary algorithms: A surve y , ” ACM Computing Surveys (CSUR) , vol. 45, no. 3, p. 35, 2013. [21] M. Gong, C. Chen, Y . Xie, and S. W ang, “Community preserving network embedding based on memetic algorithm, ” IEEE T ransactions on Emerging T opics in Computational Intelligence , no. 99, pp. 1–11, 2018. [22] W . Du, W . Y ing, P . Y ang, X. Cao, G. Y an, K. T ang, and D. W u, “Network-based heterogeneous particle swarm optimization and its application in uav communication coverage, ” IEEE T ransactions on Emer ging T opics in Computational Intelligence , pp. 1–12, 2019. [23] D. W u, N. Jiang, W . Du, K. T ang, and X. Cao, “Particle swarm optimiza- tion with moving particles on scale-free networks, ” IEEE T ransactions on Network Science and Engineering , 2018. [24] S. Das and P . N. Suganthan, “Dif ferential e volution: a survey of the state- of-the-art, ” IEEE transactions on evolutionary computation , vol. 15, no. 1, pp. 4–31, 2011. [25] Y . W ang, Z. Cai, and Q. Zhang, “Differential ev olution with composite trial vector generation strategies and control parameters, ” IEEE T rans- actions on Evolutionary Computation , vol. 15, no. 1, pp. 55–66, 2011. [26] H. Liu, Z. Cai, and Y . W ang, “Hybridizing particle swarm optimization with differential evolution for constrained numerical and engineering optimization, ” Applied Soft Computing , vol. 10, no. 2, pp. 629–640, 2010. [27] Z. Y ang, K. T ang, and X. Y ao, “Large scale evolutionary optimization using cooperative coev olution, ” Information Sciences , vol. 178, no. 15, pp. 2985–2999, 2008. [28] Z. Y ang, X. Y ao, and J. He, “Making a difference to differential ev olu- tion, ” in Advances in metaheuristics for hard optimization . Springer , 2007, pp. 397–414. [29] Y . Shang, “Global stability of disease-free equilibria in a two-group si model with feedback control, ” Nonlinear Anal Model Contr ol , vol. 20, no. 4, pp. 501–508, 2015. [30] P . V an Mieghem, J. Omic, and R. Kooij, “V irus spread in networks, ” IEEE/ACM T ransactions on Networking (TON) , vol. 17, no. 1, pp. 1–14, 2009. [31] K. Devriendt and P . V an Mieghem, “Unified mean-field framew ork for susceptible-infected-susceptible epidemics on networks, based on graph partitioning and the isoperimetric inequality , ” Physical Review E , vol. 96, no. 5, p. 052314, 2017. [32] R. Storn and K. Price, “Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces, ” Journal of global optimization , vol. 11, no. 4, pp. 341–359, 1997. [33] B.-C. W ang, H.-X. Li, J.-P . Li, and Y . W ang, “Composite differential ev olution for constrained evolutionary optimization, ” IEEE Tr ansactions on Systems, Man, and Cybernetics: Systems , 2018. [34] X. Y ao, Y . Liu, and G. Lin, “Evolutionary programming made faster , ” IEEE T ransactions on Evolutionary computation , vol. 3, no. 2, pp. 82– 102, 1999. 9 [35] T . T akahama and S. Sakai, “Constrained optimization by the ε con- strained differential ev olution with an archiv e and gradient-based mu- tation, ” in Evolutionary Computation (CEC), 2010 IEEE Congress on . IEEE, 2010, pp. 1–9. [36] B.-C. W ang, H.-X. Li, and Y . Feng, “ An improved teaching-learning- based optimization for constrained ev olutionary optimization, ” Informa- tion Sciences , vol. 456, pp. 131–144, 2018. [37] A.-L. Barab ´ asi and R. Albert, “Emergence of scaling in random net- works, ” science , vol. 286, no. 5439, pp. 509–512, 1999. [38] Y . Feng, L. Ding, Y .-H. Huang, and L. Zhang, “Epidemic spreading on weighted networks with adaptive topology based on infective informa- tion, ” Physica A: Statistical Mechanics and its Applications , vol. 463, pp. 493–502, 2016. [39] A. Khanafer, T . Bas ¸ar, and B. Gharesifard, “Stability properties of infected networks with low curing rates, ” in 2014 American Control Confer ence . IEEE, 2014, pp. 3579–3584. [40] Y . Shang, “ A lie algebra approach to susceptible-infected-susceptible epidemics, ” Electronic Journal of Differ ential Equations , vol. 2012, no. 233, pp. 1–7, 2012. [41] ——, “Lie algebraic discussion for affinity based information diffusion in social networks, ” Open Physics , vol. 15, no. 1, pp. 705–711, 2017. [42] Q. Duan, C. Shao, L. Qu, Y . Shi, and B. Niu, “When cooperative co- ev olution meets coordinate descent: Theoretically deeper understandings and practically better implementations, ” in 2019 IEEE Congr ess on Evolutionary Computation (CEC) . IEEE, 2019, pp. 721–730.