Quasi scale-free geographically embedded networks over DLA-generated aggregates

Physica A 523 (2019) 1286–1304

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Quasi scale-free geographically embedded networks over DLA-generated aggregates ∗

S. Salcedo-Sanz , L. Cuadra Department of Signal Processing and Communications, Universidad de Alcalá, Escuela Politécnica Superior, 28871 Alcalá de Henares, Madrid, Spain

highlights • • • •

We We We We

propose a meta-heuristic to construct geographically-embedded scale-free networks. consider fractal distribution of nodes from a DLA Scheme. show that it is possible to generate scale-free networks in this context. propose a Coral Reefs Optimization algorithm to construct the networks.

article

info

Article history: Received 7 November 2018 Received in revised form 24 February 2019 Available online 24 April 2019 Keywords: Geographically embedded networks Scale-free Complex network generation Diffusion limited aggregation

a b s t r a c t This paper proposes a novel model to construct quasi scale-free geographicallyembedded networks, in which the network’s nodes have been generated by means of a Diffusion Limited Aggregation (DLA) process, and the links connecting a new node i to others are generated using a link radius Ril based on the Euclidean distance. The network is constructed by means of an evolutionary-based algorithm, aiming at minimizing the error between the actual degree distribution of the network under construction and that of the scale free networks, k−γ . Several algorithms for generating scale-free geographically-embedded networks over regular Euclidean lattices can be found in the literature but, to the best of our knowledge, the general case over complex or fractal substrates had not been tackled up until now. Although well-known in very large complex networks (with a huge number of nodes and links), the scale-free property has received much less attention for small geographically-embedded networks, in which the study of networks’ properties is much more difficult. The idea of this work is to evolve the link radii for all the nodes in the network, aiming at finally fulfilling the scale-free property, if possible. Our experimental work shows that the proposed model is actually able to generate quasi scale-free geographically embedded networks in an efficient way. Discussions on the algorithm’s performance to generate scale-free networks, a special encoding to improve the search, and the algorithm’s computational evolution are given in the paper. Alternative possibilities of distribution objective (such as Poissonian, which leads to random geographically-embedded networks) are also tested and discussed in this work. © 2019 Elsevier B.V. All rights reserved.

1. Introduction What do systems as different as power grids, ecosystems, the world-wide web or a genetic network have in common? All of them can be described in terms of graphs [1]: a node represents an entity (generator/load in a power grid [2–4], ∗ Corresponding author. E-mail address: [email protected] (S. Salcedo-Sanz). https://doi.org/10.1016/j.physa.2019.04.060 0378-4371/© 2019 Elsevier B.V. All rights reserved.

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1287

different species in ecosystems [5], a web page in world-wide web [6] or a gene in a genetic network [7]), which is linked with others [8] (by electrical cables in the power grid, trophic relationships in ecosystems, protocol connections in the Internet, or gene relationships in a genetic network). These and other dissimilar systems are called complex systems, because it is extremely difficult to deduce their emerging collective behavior (for instance, synchronization [9]) from only the isolated components of the system [1]. Their topological and dynamical features can be, however, studied using the Complex Networks (CN) science [1,8,10]. The study of different complex systems using CN science has led to the identification of a number of networks types, of great interest for the description of many different systems [10–12]. Among them, ‘‘random graphs’’ have been profusely used in the description of systems such as social networks. First introduced by Erdös and Rényi [13], this type of networks are characterized by a network degree (number of links per node, k), which usually follows a Poisson distribution, though other types of distribution are also possible [14]. The second type of network models which have been deeply studied in the literature are the so-called ‘‘small-world’’ networks. These were first proposed by Watts and Strogatz in the study of social networks properties [15]. Small-world networks show simultaneously the important feature of ‘‘clustering’’, meaning that nodes that have common neighbors exhibit a higher probability of being themselves connected, and also the small-world phenomenon, meaning that their topological diameter slowly increases with the network size. Finally, Barabási and Albert introduced the concept of scale-free networks [16]. These networks are characterized by the fact that their degree distribution follows a power law P(k) ∼ k−γ . In scale-free networks, there is a large number of small degree nodes, most of which are absent in random networks, and also the probability of observing a high-degree node (hubs) is several orders of magnitude higher than in a random network [17]. Some real systems can be perfectly described as scale-free networks such as the Internet [18], metabolic networks [19], epidemic spreading [20,21], among many others [22,23]. Note that in many of these CNs (for instance, citation networks or the Internet) the position of nodes in the physical space plays no role at all [24], and the distance between two nodes is the minimum number of links that must be traveled to go from one node to the other (geodesic distance). However there are other CNs (such as transportation, infrastructure and wireless communication networks [24]) in which nodes and links are embedded in space. In this particular kind of CNs, called geographically-embedded (also spatially-embedded or simply spatial networks), nodes are located in a space associated with a metric, usually the Euclidean distance [25]. Spatial networks are usually classified into two categories: the first one, called planar networks, are those networks that can be drawn in the plane in such a way that their links do not intersect [24]. The second one involves spatial non-planar networks (for instance airline networks, cargo ship networks, or power grids) where links (which can intersect in the plane) have a cost related to their length. Although the scale-free property is well-known in very large, non-spatial complex networks (with a huge order or number of nodes), however it is not the case in small geographically-embedded networks. This is because, in spatial networks, when geometric constraints are very strong or when the cost associated with the addition of new links is large (water and gas distribution networks, power grids, or communication networks), the appearance of hubs and the scale-free feature become more difficult [24]. There are very few works dealing with the study of geographically-embedded network construction models. In [26], the authors, using as a previous layout the squares of a Sierpinski gasket (embedded in a square lattice, distorting thus its original shape), added new long-range links with a length r taken from the distribution P(r) ∼ r −α , leading to a small world nature. In [27] a model to construct random (Poissonian degree distributed) geographically-embedded networks over regular square lattices was discussed. The generating model proposed is simple, in such a way that each node is connected to others depending on a power law of the distance r between nodes, P(r) ∼ r −δ . The results obtained with this model show an efficient construction of networks, in which different regimes of interest (small-world, non-zero clustering and large world networks) were identified depending on δ . In [28] the authors proposed a model to construct also scale-free networks over regular Euclidean lattices. In this case, the embedding is driven by a constraint related to the minimization of the total length of the links in the network construction. In [29] an alternative model to generate scale-free networks over regular Euclidean lattices was proposed. This model starts from the lattice nodes and defines a ‘‘radius of action’’, R, so all the nodes within a distance R are connected to the current node. The distribution of the radii is then taken so that the number of links follows the power law distribution. This process is possible because of the regular distribution of nodes in the network. A key point to note is that all these previous approaches start from a set of nodes already distributed in a lattice, so that the links are added until the desired type of network is achieved. However, in the case of scale-free geographically-embedded networks, the aforementioned models are not applicable over random or fractal substrates (nodes layouts). Regarding this, in this paper we propose a model to generate quasi scale-free geographically-embedded networks, in which both nodes – generated by Diffusion Limited Aggregation (DLA) process [30] – and links are dynamically generated, leading to a fractal layout, unlike the papers mentioned above. Specifically, we propose a simple generation model based on assigning a link radius to each new node i, Ril , generated in the DLA aggregate. Note that this concept is similar to the radius of action defined in [29], but differs in the fact that, in [29], the nodes are previously distributed in a square lattice. A key point to note in our model for network construction establishes that each link radius may be different for each node, and it is fully related to the network generation: when a node is randomly generated by means of DLA processes, such a node is linked with all the other existing nodes in the network which are within a distance of its link radius. From now on, whenever we use the word ‘‘distance’’ in this paper, we refer to the Euclidean distance (and not the geodesic distance used in networks that are not embedded in space). In order to obtain a scale-free distribution for this growing network,

1288

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

we evolve the model as follows. We use an evolutionary-based algorithm aiming at assigning link radii to all the nodes appearing in the network, in such a way that the final degree distribution is a power law k−γ . We state this problem as an optimization task, with discrete-based encoding, in which a meta-heuristic search must be applied (since brute-force schemes are discarded due to excessive computational cost). Specifically, we consider nodes generation by a class of DLA procedure, in order to randomly generate network nodes with a fractal layout. We evaluate the performance of the Coral Reefs Optimization algorithm with Substrate Layer (CRO-SL), as evolutionary-based approach for this problem of complex networks evolution. We will show that the CRO-SL is able to generate geographically-embedded quasi scale-free networks. Different issues related to the problem encoding in the algorithm, its convergence properties and its computational cost versus search capabilities are also described in the paper. Alternative results showing that this procedure is also able to generate networks with different degree distribution such as Poisson (random network) over DLA aggregates are finally discussed. 2. Generation of geographically-embedded nodes with diffusion limited aggregation As mentioned before, one of the contributions of our proposal is based on the fact that each node is randomly generated and located (forming a fractal structure guided by a DLA), links being formed aiming at fulfilling some geometric restrictions and degree distributions (to obtain, for instance, scale free spatial network over the fractal layout). The purpose of this section is just to describe the DLA our model is based on. DLA [31] is a simple algorithm which generates random fractal clusters and aggregates [32]. An approach that hybridizes DLAs and Lindenmayer Systems (L-Systems) has been proposed very recently in [33] as a novel combinative approach to growth fractal structures. Aiming at illustrating the DLA method, we have considered the specific DLA model implementation described in [34]. This DLA model takes into account a unique seed at the beginning of the simulations, located at the center of a given discrete lattice, which will be the initial cluster to be grown. Then, particles are sequentially released at a circle distant from the cluster. The position of these particles are modified by means of a random walk. The distance between the center of the lattice and the launching circle is denoted as Ro . The random walks are defined as follows:

(

xn+1 yn+1

)

( =

xn yn

)

( +

cos(ϕ + λθn ) sin(ϕ + λθn )

) (1)

where xn and yn are the particle position at the nth step of the random walk, ϕ ∈ [−π, π] is a random angle that defines the direction of the particle’s trajectory, λ ∈ [0, 1] is a parameter that introduces the random component of the trajectories, and θn ∈ [−π, π] is a random direction. Note that λ is defined at the beginning of the random walks, whereas θn receives a random value for each walk step n. The random walk of a given particle ends if the particle visits a site neighboring the current cluster (it joins finally to this site). On the other hand, the random walk also ends if the distance between the particle and the cluster is larger than a killing radius Rk . In this latter case, the particle is eliminated and a new one is released at the launching circle Ro . Usually, the relationship between Ro and Rk can be obtained from: Ro = Rmax + Rx ,

(2)

where Rmax is the maximum distance from the center of the lattice of a particle belonging to the cluster, and Rx is a defined radius. The killing radius Rk is then defined as Rk = 2Rmax + Rx = Rmax + Ro ,

(3)

which gives the relationship between Ro and Rk . For the sake of clarity and illustrative purposes, Fig. 1 represents an example of the simulation construction by means of circles of radius Rmax , Ro and Rk . In this figure we have represented two particles, labeled ‘‘1’’ and ‘‘2’’, respectively, along with their random walk trajectories. Note that the random walk of Particle ‘‘1’’ ends successfully by reaching the cluster (particle in green aggregated to the cluster). On the contrary, Particle ‘‘2’’ has an unsuccessful random walk (particle trajectory killed after reaching circle with radius Rk ). Fig. 2 shows an example of a fractal structure generated by a DLA simulation. Note that in this case, the DLA simulation only considers integer locations for particles, so the aggregate is being constructed by direct attachment of particles to another particle already in the aggregate. Mathematically, this means that the DLA random walk equations are calculated as follows:

(

xn+1 yn+1

)

( =

xn yn

)

( +

[cos(ϕ + λθn )] [sin(ϕ + λθn )]

) (4)

where [·] stands for the nearest integer operator. However, note that this way of constructing the DLA aggregate does not lead to an adequate node set to construct a network over it. For this, Eq. (1) is preferable, and we can define an attachment radius, Ra , in such a way that a given particle is set into the DLA aggregate as soon as it finds an existing particle within a distance Ra from it. Fig. 3 shows an example of the attachment procedure to generate DLA aggregates in which the minimum distance between one particle and at least another one in the aggregate is smaller than Ra . (See Fig. 4.)

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1289

Fig. 1. Simple example to illustrate the DLA concept and the meaning of Rmax , Rl and Rk in a DLA simulation. Particles labeled ‘‘1’’ and ‘‘2’’ are released from the launching circle Rl . Particle ‘‘1’’ (green) results in being aggregated to the cluster. Particle ‘‘2’’ never reaches the cluster and ends up disappearing in Rk .

Fig. 2. Example of a fractal structure obtained by means of a DLA simulation.

Fig. 3. Example of particle attachment to the DLA aggregate considering an attachment radius Ra . See the main text for further details.

1290

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 4. Example of a DLA aggregate obtained with Eq. (1) with Ra = 10, to be used as nodes for a complex network.

3. Generating links: Growing geographically-embedded networks over DLA-distributed nodes Let us consider a model for growing networks using DLA-distributed aggregates (network nodes), in the way described in Fig. 3. The idea is to grow the network as the DLA-based nodes are being attached to the DLA aggregate. Note that this way we consider a random location for the nodes (depending on the DLA simulation), and the network is completely grown from scratch. Without loss of generality, we can consider DLA simulations starting with a single particle located at (0, 0). Then, the DLA simulation starts to add particles to the aggregate, and each time a particle is aggregated, it links to other particles already in the DLA aggregate. In order to do this, we consider an extremely simple model, based on a link radius Rl . Note that this radius is a characteristic of the particle being simulated i, so it could be better described as Ril . Note also that the number of links established if the particle i is finally attached to the current aggregate will only depend on it (Ril ). Let L be the matrix of links so that Lij stands for a binary variable describing whether or not there is a link between node i and an alternative node j. Then, each time that a particle i of the DLA is attached to the DLA aggregate, it establishes possible links to other particles already attached, in the following way:

{ Lij =

1 0

if d(i, j) < Ril otherwise

(5)

where d(i, j) stands for the distance between particle i and any other particles already on the aggregate (j). Note that if we want to ensure that all the nodes are connected with at least one other node in the network, then Ril ≥ Ra . For illustrative purposes, let us consider some examples: Fig. 5(a) shows a DLA based network with Ra = 10 and Ril = 10 (for all particles in the DLA aggregate). Note that since Ril = Ra , each particle will be attached to a very reduced number of other existing particles in the aggregate. If we keep Ra = 10 for all the DLA simulation, but Ril takes values in [10, 15, 20, 30], depending on the particle simulated, then we obtain the network in Fig. 5(b).

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1291

Fig. 5. Example of geographically embedded DLA complex networks generated with the proposed simple model considering Ra and Ril ; (a) Example of complex network with Ra = 10 and Ril = 10; (b) Example of complex network with Ra = 10 and Ril ∈ [10, 15, 20, 30].

4. Constructing quasi-scale free geographically embedded DLA-based networks The question that arises is whether or not the proposed model for network construction over geographically embedded DLA nodes can generate scale-free networks. Note that in this case the random situation of the nodes makes impossible to obtain an exact solution such as the ones proposed for square lattices in [28,29]. The approach, therefore, should be stochastic due to the nature of the considered DLA aggregates generation, and approximate solutions could arise. 4.1. Optimization problem Let us consider a DLA-based aggregate, with N nodes. This means that, after the complete DLA simulation, there will be N particles in the aggregate. Recall that the network is being constructed as the DLA simulation goes on, so each time a particle i is simulated and attached to the DLA aggregate, a new node is included in the network and the number of links in the network (L) is modified to include the new node. Let us consider a given Ra for the complete DLA simulation (all the particles will be attached to the DLA aggregate with the same Ra ) and specific Ril radius for each particle (Ra ≥ Ril ). Let r = [R1l , . . . , RNl ] the radii associated with the N particles forming the DLA aggregate. The idea is to obtain a vector r∗ which makes the network degree to have a scale-free distribution, i.e., such that it minimizes the following objective function: f (r) =

kmax ∑ ⏐ ⏐ ⏐pk (r) − k−γ ⏐

(6)

k=2

where pk stands for the degree distribution of the DLA-based network obtained with a vector r, and kmax is the maximum degree considered depending on the network size. That is, the objective function f (r) quantifies to what extent the

1292

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 6. Compressed encoding example (β = 5), useful in the evolution of complex networks. Table 1 Parameters of the CRO-SL used in the evolution of the networks considered. See [37,38] for further details about the parameters. CRO-SL

Parameters

Initialization External sexual reproduction Substrates Internal sexual reproduction Larvae setting Asexual reproduction Depredation Stop criterion

Reef size = 50 × 40, ρ0 = 0.9 Fb = 0.80 T = 5 substrates: HS, DE, 2Px, MPx, SM 1 − Fb = 0.20 κ=3 Fa = 0.05 Fd = 0.15, Pd = 0.05 tmax = 500 iterations

distribution function of the under-construction spatial network, pk (r), differs from that of a scale-free one, k−γ . Another point to note is that we aim to find out whether or not there exists a r∗ which leads to a power law distribution in the network with a given γ . This problem is therefore stated as an integer optimization problem, in which the final network degree distribution will completely depend on r. The problem is discrete, highly non-linear, and the size of the search space is huge when N grows, which discards exact solutions via brute force algorithms. In these kind of problems, meta-heuristics approaches such as Evolutionary Computation-based algorithms are able to obtain very good solutions with a moderate computational complexity. We therefore propose to apply a kind of Evolutionary Algorithm to solve this optimization associated with Scale-Free DLA-based networks. 4.2. Evolutionary algorithm proposed Evolutionary algorithms are meta-heuristics approaches, mainly focused on solving optimization problems, based on mimicking different process of natural evolution in a computer [35,36]. In an evolutionary-type algorithm, a population of tentative solutions to the problem is maintained in the algorithm, usually encoded as vectors of several variables, and different operators are sequentially applied to it, in order to obtain better solutions to the problem. This process is guided by an objective function, which must be optimized (minimized or maximized) at the end of the algorithm. In the most basic version of any evolutionary algorithm, the evolution operators are usually a kind of crossover and mutation to generate new solutions (algorithm’s exploration) and some kind of selection for promoting the best solutions found so far in the evolution and discarding the worst ones. Depending on the problem’s encoding, these operators may suffer modifications, or even need complicated versions to ensure feasibility of solutions, etc. 4.2.1. The proposed coral reefs optimization with substrate layer The Coral Reef Optimization algorithm (CRO) [37] (further described in [38]), is an evolutionary-type algorithm based on the behavior of the processes occurring in a coral reef. Additionally, in [39], a new version of the CRO algorithm with Substrate Layer (CRO-SL) has been presented. In the CRO-SL approach, several substrate layers (specific parts of the population) have been introduced in such a way that each substrate may represent different processes (different models, operators, parameters, constraints, repairing functions, etc.) useful in the evolution. Specifically, in [40] a version of the CRO-SL algorithm has proposed, in which each substrate layer represents a different search procedure, leading to a co-evolution competitive algorithm. This version of the CRO-SL has been successfully tested in different applications and problems such as micro-grid design [41], vibration cancellation in buildings [42,43], or in the evaluation of novel non-linear search procedures [44]. This is also the CRO-SL algorithm version used in this paper. Regarding the algorithm’s encoding for the optimization problem at hand, we consider integer vectors as solutions, x ∈ N. Note that using this encoding, the length of each individual is l = N. This encoding provides a compact version of the algorithm, and allows using some different searching procedures such as Harmony Search or Differential Evolution. The main problem with a direct encoding of N integer values in the CRO-SL algorithm is that, as N grows, the searching capabilities of the algorithm can be affected, since the search space is huge. It is possible to manage shorter encodings

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1293

Fig. 7. DLA-generated aggregates, i.e. (fractal layouts of N nodes spatially embedded in the space) with Ra = 10; (a) N = 100; (b) N = 200; (c) N = 400; (d) N = 800.

by using a compressed version of the encoding, in such a way that each element of the encoding represents β actual values, such as we proposed in [45]. Fig. 6 shows an example of this compressed encoding, which reduces the current encoding length l = N to l′ = Nβ . Of course, the resolution of the search space is smaller than in the original encoding when the compressed encoding is applied, but on the other hand, it is expected that the CRO-SL algorithm searches for better solutions in this smaller search space. The considered substrates for solving the stated problem are detailed below. Note that they are general purpose substrates, such as Differential Evolution or Harmony Search-based, and other specific substrates with crossovers and mutations adapted to the chosen encoding. Five different substrates will be described and evaluated later in the experimental section.

• Differential Evolution-based operator (DE): This operator is based on the evolutionary algorithm with that

•

•

•

•

name [46], a method with powerful global search capabilities. DE introduces a differential mechanism for exploring the search space. Hence, new larvae are generated by perturbing the population members using vector differences of individuals. Perturbations are introduced by applying the rule x′i = x1i + F (x2i − x3i ) for each encoded parameter on a random basis, where x′ corresponds to the output larva, xt are the considered parents (chosen uniformly among the population), and F determines the evolution factor weighting the perturbation amplitude. Harmony Search-based operator (HS): Harmony Search [47] is a population based meta-heuristic that mimics the improvisation of a music orchestra while its composing a melody. This method integrates concepts such as harmony aesthetics or note pitch as an analogy for the optimization process, resulting in a good exploratory algorithm. HS controls how new larvae are generated in one of the following ways: (i) with a probability HMCR ∈ [0, 1] (Harmony Memory Considering Rate), the value of a component of the new larva is drawn uniformly from the same values of the component in the other corals. (ii) with a probability PAR ∈ [0, 1] (Pitch Adjusting Rate), subtle adjustments are applied to the values of the current larva, replaced with any of its neighboring values (upper or lower, with equal probability). Two points crossover (2Px): 2PX [35] is considered one of the standard recombination operators in evolutionary algorithms. In the standard version of the operator, two parents from the reef population are provided as input. A recombination operation from two larvae is carried out by randomly choosing two crossover points, interchanging then each part of the corals between those points. Multi-points crossover (MPx): Similar to the 2PX, but in this case the recombination between the parents is carried out considering a high number of crossover points (M), and a binary template which indicates whether each part of one parent is interchanged with the corresponding of the other parent. Standard integer Mutation (SM): This operator consists of a standard mutation in integer-based encodings. It consists of mutating an element of a coral with another valid value (different from the previous one). Note that the SM operator links a given coral (possible solution) to a neighborhood of solutions which can be reached by means of a single change in an element of the coral.

1294

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304 Table 2 Main parameters of the experiments carried out: DLA aggregate size (or number of nodes, N), the maximum degree to evaluate the scale-free property in each network (kmax ), the compression rate in the CRO-SL encoding (β ) and the final encoding length of the solutions in the CRO-SL (l′ ). N

kmax

β

l′

100 200 400 800

10 12 15 20

5 5 10 10

20 40 40 80

Table 3 Results in terms of f (r∗ ) obtained with the CRO-SL algorithm for the different networks considered. r∗ labels the vector that minimized the objective function f (r) stated by Eq. (6). N

f (r∗ ) (Eq. (6))

100 200 400 800

0.0136 0.0111 0.0140 0.0142

5. Experiments and results In the experiments that are described below, we have explored the feasibility of both the proposed model and the CRO-SL algorithm to tackle it in spatial network, with DLA aggregates sizes in the order of hundreds nodes. Spatial networks with such relatively small order (of the order of hundreds of nodes) can appear in complex problems in which the distribution of the nodes is not based on a lattice, as, for instance, in regional science analysis [48,49] or in rainfall modeling [50]. 5.1. Experimental setting In this section we show different computational results obtained with the CRO-SL in the evolution of different DLAbased networks, over DLA aggregates with N = 100, 200, 400 and 800 nodes, respectively. The resulting DLA aggregates considered have been represented in Fig. 7. In all cases, a common Ra = 10 value has been considered to form the aggregates. The design of the networks over the DLA aggregates is carried out by the CRO-SL, as described in the previous section. Note that the link radius to be assigned by the algorithm to all nodes in the DLA aggregates has been forced to fulfill the property Ra ≤ Ril ≤ 100. Table 1 shows the corresponding values for the CRO-SL parameters considered in the experiments carried out. It is very important to note that, when we aim to dynamically generate a quasi scale-free spatial network with a final number of nodes (N) distributed in space following a fractal structure from the DLA, the CRO-SL algorithm evaluates about 100,000 possible networks (500 iterations × Reef size). This value has been kept common for all the DLA aggregates considered, independently of their size. A summary of some important parameters for the experiments −including the DLA aggregate size (or number of nodes, N), the maximum degree to evaluate the scale-free property in each network (kmax ), the compression rate in the CRO-SL encoding (β ) and the final encoding length of the solutions in the algorithm (l′ )− has been listed in Table 2. 5.2. Testing the proposal for different network orders, N: Results Table 3 will assist us in discussing a first approximation to the results. It represents, as a function of the number of nodes N, the minimum value f (r∗ ) of the objective function (Eq. (6)) that the proposed CRO-SL algorithm has found for each DLA aggregate considered in the experimental setting. Taking into account that f (r) minimizes the error between both distributions (that of the network under construction, pk (r), and that of the scale-free, k−γ ), Table 3 shows that our method is able to generate geographically embedded networks over DLA aggregates with a degree distribution very close to be scale-free in all cases considered. With this in mind, we will discuss in the following paragraphs the networks obtained for N = 100, 200, 400 and 800 nodes, respectively. First, Fig. 8 shows the results obtained for the smallest DLA aggregate considered, with N = 100 nodes. It represents the node degree distribution pk (r∗ ) corresponding to the spatial network generating-method driven by the CRO-SL algorithm. This figure, like all those that follow, contains two insets. The first one, on the upper-right part of the figure, shows the resulting quasi-scale free geographically embedded network. The second inset, on the lower left part of the figure, represents, for each node, the link radius Rl of any node corresponding to the best solution obtained with the CRO-SL

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1295

Fig. 8. Example of a quasi scale-free geographically-embedded complex network with N = 100 nodes, evolved with the CRO-SL algorithm. In the node degree distribution for the network, blue circles stand for the power law distribution k−1.52 , and red points for the actual degree distribution of the obtained network. The resulting final spatial network obtained and the best solution obtained with the CRO-SL are also depicted in this figure.

Fig. 9. Example of a quasi scale-free geographically-embedded complex network with N = 200 nodes, evolved with the CRO-SL algorithm. In the node degree distribution for the network blue circles stand for the power law distribution k−1.55 , and red points for the actual degree distribution of the obtained network. The resulting final spatial network obtained and the best solution obtained with the CRO-SL are also depicted in this figure.

(r∗ ). Note that, in this network, as listed in Table 2, kmax = 10, i.e., we consider a linear bidding k ranging from 2 to 10 in Eq. (6), plus another upper value of k averaging the contribution of higher degrees. We have followed this methodology for all the networks explored in the work. The optimization process applied to this network has led to a quasi-scale-free network with a power law distribution with k−1.52 . It is interesting to check that there are a small number of nodes that act as hubs and a much larger number of nodes with a reduced number of links. This behavior is consistent in all the DLA-based networks evolved in this work. Fig. 9 illustrates the results obtained by the CRO-SL in the network evolution for the N = 200 DLA aggregate. The insets have the same meaning that those of the previous one. Fig. 9 shows the network obtained after the optimization process, which has obtained an excellent approximation of the network distribution node degree, leading to a power law distribution with k−1.55 . Note that, as shown in Table 2, we have explored a linear bidding of 12 values of the node degree k in the network, ranging from 2 to 12, plus another point in higher values of k. The geographically-embedded network corresponding to the best solution r∗ found by the CRO-SL algorithm has been represented on the bottom-right inset. Note the runs of β = 5 (Table 2) is the same as that of the N = 100 case. Fig. 10 shows the results obtained by the CRO-SL in the evolution of the network over the DLA aggregated with N = 400 nodes. In this case, we have considered a compressed encoding with β = 10, which leads to a l′ = 400 = 40, similar to 10 the N = 200 case, and larger than the N = 100 network encoding (see Table 2). We have found that this compressed encoding provides the best results. This figure shows the resulting network generated by the CRO-SL algorithm, which is constructed to very approximately follow a power law distribution k−1.59 . In this case we have explored a linear bidding with 15 values of the degree k, from 2 to 15, and a rest in higher values of k. The best solution r∗ found by the CRO-SL algorithm has also been displayed in Fig. 10. Note the runs of β = 10 equal values in the solution. For illustrative purposes, Fig. 11 shows 6 steps of the network evolution for network illustrated by the uppermost inset of Fig. 10. In Fig. 11, it is possible to see the process of network construction as the nodes are being attached to the aggregate in the DLA process. It is important to take into account that the spatial network construction depends on the

1296

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 10. Example of a quasi scale-free geographically-embedded complex network with N = 400 nodes, evolved with the CRO-SL algorithm. In the node degree distribution for the network blue circles stand for the power law distribution k−1.59 , and red points for the actual degree distribution of the obtained network. The resulting final spatial network obtained and the best solution obtained with the CRO-SL are also depicted in this figure.

Fig. 11. Network evolution process in 6 steps for the N = 400 quasi scale-free case (best solution obtained with the CRO-SL algorithm).

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1297

Fig. 12. Example of a quasi scale-free geographically-embedded complex network with N = 800 nodes, evolved with the CRO-SL algorithm. In the node degree distribution for the network blue circles stand for the power law distribution k−1.62 , and red points for the actual degree distribution of the obtained network. The resulting final spatial network obtained and the best solution obtained with the CRO-SL are also depicted in this figure.

Fig. 13. Geographically-embedded network evolution process in 5 steps for the N = 800 quasi scale-free case (best solution obtained with the CRO-SL algorithm).

position of the randomly generated nodes (we focus on geographically-embedded networks), controlled by Ra and also in the values of Ril , which are evolved by the CRO-SL algorithm. Finally, Fig. 12 shows the results obtained by the CRO-SL in the evolution of the network over the largest DLA aggregate considered with N = 800 nodes. As in the other cases, this figure shows, respectively, the degree distribution of the evolved network (quasi scale-free), which very approximately follows a power law distribution k−1.62 , the final network

1298

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 14. Network duplication (cloning and shifting) process to 2N for the network defined over the DLA aggregate with N = 200; (a) Final network; (b) Quasi scale-free distribution of the final duplicated network.

obtained, and the best solution evolved by the CRO-SL algorithm. Note that in this network kmax = 20 (see Table 2), i.e. we consider a linear biding with a range of k values from 2 to 20 in Eq. (6), plus another upper value of k averaging the contribution of higher degrees. The construction process represented in Fig. 13 shows 5 steps of the network evolution, for the best solution found by the CRO-SL. As in the previous case for N = 400, in this figure it is possible to see the process of network construction as the nodes are being attached in the DLA process. As can be seen in the results obtained, we have shown that it is possible to generate quasi-scale-free geographicallyembedded random networks over DLA aggregates, by considering a very simple model of distances between nodes. It is necessary to solve an optimization problem, which is hard, since it must optimize the link radius of all the randomly generated nodes which form the network. We have shown how the CRO-SL algorithm with compressed encoding is able to successfully solve this task, finding near optimal solutions to the optimization problem. 5.3. Spatial duplication of quasi scale-free networks. results Since we are dealing with geographically embedded networks, it may have sense to generate larger networks starting from an initial one, by means of spatial duplication (i.e. cloning the initial network) and shifting to a different position, while linking it somehow to the original one. The result should be another network with the double of nodes. This process can be iterated in order to obtain even larger networks. The question arising here is whether or not these structures can preserve the quasi scale-free property, which will depend on how each cloning network links to the original one. In order to illustrate this process, Fig. 14 shows the spatial duplication (cloning and shifting to the right) of the N = 200 network previously obtained with the CRO-SL (see Fig. 9), and its final degree distribution. The link between the cloned network and the initial one has been assigned in this case by brute force, in such a way that Eq. (6) is minimized. In cases with larger number of cloning networks this problem is itself NP-hard, so we have obtained the final links between networks by inspection in a reduced number of tries. As a first example, we have duplicated and shifted the spatial network with N = 200 nodes, leading to the spatial network represented in Fig. 14(a), with N = 2 × 200 = 400 nodes. As shown in Fig. 14(b), the resulting network keeps very approximately the scale free property.

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1299

Table 4 Results in terms of f (r∗ ) (Eq. (6)) obtained with the CRO-SL algorithm for the different networks considered, with D = N, 2N, 4N and 8N. N

f (r∗ )

γ

2N

f (r∗ )

γ

4N

f (r∗ )

γ

8N

f (r∗ )

γ

100 200 400 800

0.0136 0.0111 0.0140 0.0142

1.51 1.55 1.59 1.62

200 400 800 1600

0.0132 0.0135 0.0132 0.0137

1.52 1.55 1.59 1.62

400 800 1600 3200

0.0136 0.0137 0.0128 0.0140

1.51 1.55 1.59 1.62

800 1600 3200 6400

0.0160 0.0140 0.0132 0.0148

1.52 1.56 1.59 1.62

Fig. 15. Duplicated networks (4N) obtained for all the DLA aggregates considered; (a) N = 100; (b) N = 200; (c) N = 400; (d) N = 800.

We can analyze this duplication process in detail by using the networks obtained in all the DLA aggregates considered in this paper. We have obtained duplicated networks (cloning and shifting to the right) up to 8 × N nodes for all the networks obtained in the previous section. Table 4 shows the results obtained for all the DLA-based networks in the duplication process. The results obtained show that the duplication process does not affect much to the scale-free property, as expected. The value of γ is the same or very similar for all the duplicated networks, and we can obtain even better networks in terms of Eq. (6) than the initial one in some cases. For example, in the case in which the network to be duplicated is the one with N = 100 nodes, the duplication to 2 × N results in a network with a value of Eq. (6) of 0.0132, better than the initial 0.0136 obtained with the CRO-SL evolution. We have can see some other cases, such as the duplication of the N = 400 DLA network to 2 × N, 4 × N and 8 × N, or the duplication of the N = 800 network to 2 × N and 4 × N. We could not find a better network by duplication in the case of N = 200, since the initial network was of very high quality. In this case, network duplication leads to larger networks with worse value of the objective function. Anyway, it is convenient to recall that the initial networks are obtained by a meta-heuristic algorithm, which does not guarantee optimality of the result, but a good-enough behavior, which may explain the improvement of the network quality by spatial duplication. Fig. 15 shows the 4 × N duplication of all the DLA networks generated in this paper. As can be seen in Fig. 16, all of them keep following a quasi scale-free degree distribution. 5.4. Generation of DLA-based quasi random spatial networks. results We finally evaluate the feasibility of the proposed approach to evolve DLA-based networks following a Poisson distribution (random networks) instead of scale-free networks. For this, we just change Equation (6) by the following

1300

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 16. Degree distribution of the 4 × N duplication process for all the DLA aggregates considered: (a) N = 100; (b) N = 200; (c) N = 400; (d) N = 800.

objective function: f (r) =

⏐ kmax ⏐ (−λ) k⏐ ∑ ⏐ ⏐pk (r) − e (λ) ⏐ ⏐ ⏐ k!

(7)

k=2

where λ is a parameter of the Poisson distribution to be obtained for each case. Fig. 17 shows the best quasi-random network obtained with the CRO-SL algorithm for the case of the DLA aggregate with N = 400 nodes. As shown, the degree distribution in this case approaches a Poisson distribution with λ = 9. The network obtained is quasi-random, without any kind of hubs. The differences between the two r solutions obtained for the quasi scale-free (Fig. 10) and the quasi-random network (Fig. 17) are significant. Fig. 18 shows the final network construction for this case of quasi-random N = 400 network, in six steps. Fig. 18 help us compare this process with that corresponding to the quasi-scale free solution (Fig. 11), and check out the differences in the network construction from the beginning of the process. 5.5. Computational complexity and algorithm’s convergence: Discussion We have shown that the proposed model based on evolutionary computing is able to generate quasi scale-free geographically embedded networks over DLA aggregates in an efficient way. The evolved networks are of relative small size (from 100 to 800 nodes), though we can obtain larger networks applying spacial duplication and shifts over the initial quasi scale-free networks. The reason of this small order of the networks explored is the computational cost of the meta-heuristics, which must encode link radii for all the nodes in the network. We have applied a compressed encoding in the CRO-SL algorithm aiming to alleviate the computational burden of the approach and to obtain a search of higher quality. Note that the compressed encoding with factor β implies that groups of β nodes will have the same link radii assigned (see solutions obtained by the CRO-SL algorithm in Figs. 8–12. The improvement in performance of the algorithm when the compressed encoding is included indicates that the initial search space is huge, which causes problems even

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1301

Fig. 17. Example of a quasi random (Poissonian) geographically-embedded network with N = 400 nodes, evolved with the CRO-SL algorithm. In the (−9) k

node degree distribution for the network, blue circles stand for the distribution e k! 9 , while red points represent those of actual degree distribution of the generated network. The resulting final spatial network obtained and the best solution obtained with the CRO-SL are also depicted in this figure.

Fig. 18. Spatial network evolution process in 6 steps for the N = 400 quasi-random network case (best solution obtained with the CRO-SL algorithm).

for small order networks cases. Regarding this, Fig. 19 shows the CRO-SL convergence process (fitness function evolution over 500 generations), for all the networks considered. All the sub-figures represent the best solution in each generation of the algorithm. Note that there are small differences in the convergence process depending on the network. Though the differences between these figures are small, it seems that the algorithm converges faster for the case N = 100, whereas

1302

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

Fig. 19. Convergence process of the CRO-SL fitness function as a function of the number of generations, for all the networks considered, with order: (a) N = 100; (b) N = 200; (c) N = 400; (d) N = 800.

in the N = 800 network the convergence is slower. Recall that in all cases, however, the results in terms of the objective function is very good, providing a quasi scale-free network. Alternatively to the compressed encoding used in this paper, there are other possibilities not explored in this paper, but open for future work. For example, we could use different patterns of compression instead of repeating the same element β times. In other words, we could substitute the current compressed encoding 20 → 20 20 20 20 20 (β = 5), for another pattern of compression such as 20 → 20 19 18 22 19 (β = 5, with elements 2 values over and under the central value (20)). Of course very different compression patterns could be defined. We could also construct compressed encodings using L-Systems [51], which could be in turn optimized with a meta-heuristic approach. To complete this computational and convergence analysis, Fig. 20 shows the larvae evolution in the CRO-SL for the problem over the DLA aggregate with N = 400 (evolution towards a scale-free network). Fig. 20(a) represents the percentage of times that a given substrate produces the best larva (solution) in a generation. In turn, Fig. 20(b) shows the number of larvae which are located into the reef in each generation. These two figures provide a good idea of the behavior of each substrate in the CRO-SL algorithm. As can be seen in Fig. 20(a), the 2Px and MPx operators seem to be the most effective in producing the best solutions over the generations. Standard mutation (SM) and DE seem to help less in the algorithm’s search, but have some contribution. The HS operator seems to be the less operational one in this problem. Fig. 20(b) provides a more reliable information on the algorithm’s evolution, since it considers the larvae (solutions) which are actually attached to the reef. This figure makes clear that the MPx substrate is the one that help the best algorithm’s in the searching process. In the first stages of the search, all operators seems to contribute to the algorithm’s search, but at the final part of the evolution, the contribution of HS and DE operators is small. Fig. 20(b) also shows that the algorithmic behavior of the CRO-SL is the expected for a well-designed meta-heuristic, since the number of new solutions attached to the reef is very large at the beginning, and small at the end of the evolution, when the algorithm is almost converged, leading to a good exploration of the search space. Finally, note that small size geographically-embedded networks (with hundreds or few thousands nodes) appear in different real problems such as regional science analysis [48,49] or rainfall modeling [50], among others. Thus, proposing

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1303

Fig. 20. Percentage of best larvae obtained from each substrate and number of new larvae which are able to get into the reef per generation and substrate, for the DLA aggregate evolution of N = 400 to a scale-free network; (a) Percentage of best larvae formed; (b) Number of new larvae in the reef.

models which show that these networks, defined over non-trivial substrates such as DLA aggregates, can be quasi scale-free is of importance in the application of CN to these and other similar problems. 6. Conclusions In this paper we have proposed both a novel model to construct quasi scale-free geographically-embedded networks and the fully tailored evolutionary algorithm to tackle its complexity. On the one hand, the proposed geographicallyembedded network generating model exhibits the novelty of adding both nodes (generated by a Diffusion Limited Aggregation (DLA) process) and links (by connecting each new node i to others within a link radius Ril based on the Euclidean distance). Nodes and links are generated aiming at minimizing the error between the degree distribution of the network under construction and that of the scale free networks, k−γ . To our knowledge, this is the first time that general geographically-embedded networks over complex substrates, different from those in lattices, have been considered for scale-free synthesis. On the other hand, since this generation problem has been stated as a minimization task in which brute-force schemes have been discarded (because of their excessive computational cost), we have used an evolutionarybased algorithm, a Coral Reefs Optimization algorithm with Substrate Layer (CRO-SL). This evolutionary algorithm works by assigning a link radius to all the nodes in the network, in such a way that the final degree distribution is (quasi) scale-free. This optimization task has been shown to lead to a large encoding of the problem, which can be alleviated by introducing a compressed encoding factor. This improves the searching capabilities of the CRO-SL algorithm. We have tested the proposal in several DLA-based aggregates with 100, 200, 400 and 800 nodes, obtaining quasi-scale free networks in all cases. This shows that it is possible to define geographically-embedded networks with the scale-free property even in complex/fractal scenarios such as nodes coming from DLA aggregates. We have also shown that the proposed approach is flexible enough to generate networks with alternative degree distributions, such as, for instance, the Poissonian one of random networks, by simply changing the objective function of the algorithm. This work opens up new research possibilities and lines, in different fields. As previously mentioned, the compressed encoding is a promising technique to tackle large-scale optimization problem, but can be improved, by considering

1304

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

different compression patterns or even L-Systems. Another interesting research line is related to the possibility of defining other types of models for network construction, such as ‘‘competition’’ networks, useful for modeling competition among different companies or firms in a zone. These competition network models can be defined by considering not only the link radius of the appearing node in the DLA, but also the link radii of already existing nodes in the aggregate. In terms of applications, Renewable Energy or Regional Science problems are quite appealing. For instance, studying the relationships among wind farms in a country using spatial complex networks could be extremely interesting. Also, there are lots of problems in Regional Science related to transport networks, economic relationships among cities, etc. Acknowledgment This work has been partially supported by Ministerio de Economía y Competitividad of Spain (Grant Ref. TIN201785887-C2-2-P). References [1] A.-L. Barabási, Network Science, Cambridge University Press, 2016. [2] W. Guo, H. Wang, Z. Wu, Robustness analysis of complex networks with power decentralization strategy via flow-sensitive centrality against cascading failures, Physica A 494 (2018) 186–199. [3] C. Chen, X. Zhou, Z. Li, Z. He, Z. Li, X. Lin, Novel complex network model and its application in identifying critical components of power grid, Physica A 512 (2018) 316–329. [4] L. Cuadra, M.D. Pino, J.C. Nieto-Borge, S. Salcedo-Sanz, Optimizing the structure of distribution smart grids with renewable generation against abnormal conditions: A complex networks approach with evolutionary algorithms, Energies 10 (8) (2017) 1097. [5] B.D. Fath, W. Grant, Ecosystems as evolutionary complex systems: Network analysis of fitness models, Environ. Model. Softw. 22 (5) (2007) 693–700. [6] A.-L. Barabási, R. Albert, H. Jeong, Scale-free characteristics of random networks: the topology of the world-wide web, Physica A 281 (1–4) (2000) 69–77. [7] C. Christensen, A. Gupta, C.D. Maranas, R. Albert, Large-scale inference and graph-theoretical analysis of gene-regulatory networks in B. Subtilis, Physica A 373 (2007) 796–810. [8] L. Lü, T. Zhou, Link prediction in complex networks: A survey, Physica A 390 (6) (2011) 1150–1170. [9] A. Arenas, A. Díaz-Guilera, J. Kurths, Y. Moreno, C. Zhou, Synchronization in complex networks, Phys. Rep. 469 (3) (2008) 93–153. [10] M.E. Newman, The structure and function of complex networks, SIAM Rev. 45 (2) (2003) 167–256. [11] S.H. Strogatz, Exploring complex networks, Nature 410 (6825) (2001) 268. [12] R. Albert, A.-L. Barabási, Statistical mechanics of complex networks, Rev. Modern Phys. 74 (1) (2002) 47. [13] P. Erdos, A. Rényi, On random graphs i, Publ. Math. Debrecen 6 (1959) 290–297. [14] M.E. Newman, S.H. Strogatz, D.J. Watts, Random graphs with arbitrary degree distributions and their applications, Phys. Rev. E 64 (2) (2001) 026118. [15] D.J. Watts, S.H. Strogatz, Collective dynamics of ‘‘small-world’’ networks, Nature 393 (6684) (1998) 440. [16] A.-L. Barabási, R. Albert, Emergence of scaling in random networks, Science 286 (5439) (1999) 509–512. [17] A.-L. Barabási, E. Bonabeau, Scale-free networks, Sci. Am. 288 (5) (2003) 60–69. [18] R. Albert, H. Jeong, A.-L. Barabási, Internet: Diameter of the world-wide web, Nature 401 (6749) (1999) 130. [19] H. Jeong, B. Tombor, R. Albert, Z.N. Oltvai, A.-L. Barabási, The large-scale organization of metabolic networks, Nature 407 (6804) (2000) 651. [20] J. Liu, T. Zhang, Epidemic spreading of an SEIRS model in scale-free networks, Commun. Nonlinear Sci. Numer. Simul. 16 (8) (2011) 3375–3384. [21] T. Li, Y. Wang, Z.-H. Guan, Spreading dynamics of a SIQRS epidemic model on scale-free networks, Commun. Nonlinear Sci. Numer. Simul. 19 (3) (2014) 686–692. [22] A.-L. Barabasi, Linked: How Everything is Connected to Everything Else and What It Means, Plume. Basic Books, 2003. [23] R. Pastor-Satorras, A. Vespignani, Epidemic dynamics and endemic states in complex networks, Phys. Rev. E 63 (6) (2001) 066117. [24] M. Barthélemy, Spatial networks, Phys. Rep. 499 (1–3) (2011) 1–101. [25] M. Barthelemy, Morphogenesis of Spatial Networks, Springer, 2018. [26] M.R. Roberson, D. Ben-Avraham, Kleinberg navigation in fractal small-world networks, Phys. Rev. E 74 (1) (2006) 017101. [27] K. Kosmidis, S. Havlin, A. Bunde, Structural properties of spatially embedded networks, Europhys. Lett. 82 (4) (2008) 48005. [28] A.F. Rozenfeld, R. Cohen, D. Ben-Avraham, S. Havlin, Scale-free networks on lattices, Phys. Rev. Lett. 89 (21) (2002) 218701. [29] C.P. Warren, L.M. Sander, I.M. Sokolov, Geography in a scale-free network model, Phys. Rev. E 66 (5) (2002) 056105. [30] H. Stanley, A. Bunde, S. Havlin, J. Lee, E. Roman, S. Schwarzer, Dynamic mechanisms of disorderly growth: Recent approaches to understanding diffusion limited aggregation, Physica A 168 (1) (1990) 23–48. [31] T.A. Witten, L.M. Sander, Diffusion-limited aggregation, Phys. Rev. B 27 (9) (1983) 5686. [32] H. Stanley, A. Coniglio, S. Havlin, J. Lee, S. Schwarzer, M. Wolf, Diffusion limited aggregation: a paradigm of disorderly cluster growth, Physica A 205 (1–3) (1994) 254–271. [33] S. Salcedo-Sanz, L. Cuadra, Hybrid L-systems–Diffusion Limited Aggregation schemes, Physica A 514 (15) (2019) 693–700. [34] S. Ferreira Jr, S. Alves, A. Brito, J. Moreira, Morphological transition between diffusion-limited and ballistic aggregation growth patterns, Phys. Rev. E 71 (5) (2005) 051402. [35] A.E. Eiben, J.E. Smith, Introduction to Evolutionary Computing, 53, Springer-Verlag, Natural Computing Series, 2003. [36] A.E. Eiben, J. Smith, From evolutionary computation to the evolution of things, Nature 521 (7553) (2015) 476. [37] S. Salcedo-Sanz, J. Del Ser, I. Landa-Torres, S. Gil-López, J. Portilla-Figueras, The coral reefs optimization algorithm: a novel metaheuristic for efficiently solving optimization problems, Sci. World J. 2014 (2014). [38] S. Salcedo-Sanz, A review on the coral reefs optimization algorithm: new development lines and current applications, Prog. Artif. Intell. 6 (1) (2017) 1–15. [39] S. Salcedo-Sanz, J. Muñoz-Bulnes, M.J. Vermeij, New coral reefs-based approaches for the model type selection problem: a novel method to predict a nation’s future energy demand, Int. J. Bio-insp. Comput. 10 (3) (2017) 145–158. [40] S. Salcedo-Sanz, C. Camacho-Gómez, D. Molina, F. Herrera, A coral reefs optimization algorithm with substrate layers and local search for large scale global optimization, in: 2016 IEEE Congress on Evolutionary Computation (CEC), IEEE, 2016, pp. 3574–3581. [41] S. Salcedo-Sanz, C. Camacho-Gómez, R. Mallol-Poyato, S. Jiménez-Fernández, J. Del Ser, A novel coral reefs optimization algorithm with substrate layers for optimal battery scheduling optimization in micro-grids, Soft Comput. 20 (11) (2016) 4287–4300.

S. Salcedo-Sanz and L. Cuadra / Physica A 523 (2019) 1286–1304

1305

[42] S. Salcedo-Sanz, C. Camacho-Gómez, A. Magdaleno, E. Pereira, A. Lorenzana, Structures vibration control via tuned mass dampers using a co-evolution coral reefs optimization algorithm, J. Sound Vib. 393 (2017) 62–75. [43] C. Camacho-Gómez, X. Wang, E. Pereira, I. Díaz, S. Salcedo-Sanz, Active vibration control design using the coral reefs optimization with substrate layer algorithm, Eng. Struct. 157 (2018) 14–26. [44] S. Salcedo-Sanz, Modern meta-heuristics based on nonlinear physics processes: A review of models and design procedures, Phys. Rep. 655 (2016) 1–70. [45] S. Salcedo-Sanz, A. Gallardo-Antolín, J.M. Leiva-Murillo, C. Bousono-Calzón, Offline speaker segmentation using genetic algorithms and mutual information, IEEE Trans. Evol. Comput. 10 (2) (2006) 175–186. [46] R. Storn, R. storn and k. price, j. global optim. 11, 341 (1997)., J. Global Optim. 11 (1997) 341. [47] Z.W. Geem, J.H. Kim, G.V. Loganathan, A new heuristic optimization algorithm: harmony search, Simulation 76 (2) (2001) 60–68. [48] R. Mansilla, R. Mendozas, The Network of Mexican Cities, arXiv preprint arXiv:1006.3887. [49] Z. Neal, Is the urban world small? the evidence for small world structure in urban networks, Netw. Spat. Econ. (2018) 1–17. [50] S.K. Jha, B. Sivakumar, Complex networks for rainfall modeling: spatial connections, temporal scale, and network size, J. Hydrol. 554 (2017) 482–489. [51] S. Salcedo-Sanz, L. Cuadra, Hybrid L-systems-Diffusion Limited Aggregation schemes, Physica A 514 (2019) 592–605.