A parallel adaptive quantum genetic algorithm for the controllability of arbitrary networks

In this paper, we propose a novel algorithm—parallel adaptive quantum genetic algorithm—which can rapidly determine the minimum control nodes of arbitrary networks with both control nodes and state nodes. The corresponding network can be fully controlled with the obtained control scheme. We transformed the network controllability issue into a combinational optimization problem based on the Popov-Belevitch-Hautus rank condition. A set of canonical networks and a list of real-world networks were experimented. Comparison results demonstrated that the algorithm was more ideal to optimize the controllability of networks, especially those larger-size networks. We demonstrated subsequently that there were links between the optimal control nodes and some network statistical characteristics. The proposed algorithm provides an effective approach to improve the controllability optimization of large networks or even extra-large networks with hundreds of thousands nodes.


Introduction
The real world consists of ubiquitous intricate and intertwined networks. Some are tangible, such as traffic networks [1,2], power networks [3], and financial networks [4], whereas others are invisible networks that penetrate into every aspect of our lives, such as interpersonal relationship networks [5,6], wireless networks [7,8], and ecological networks [9]. The expected goal of research into a complex network is to be able to regulate and control it from the outside, achieve the desirable state or performance by injecting outward control signals to some network nodes (called driver nodes), and ultimately achieve the real controllability of a complex network.
Many studies have been conducted on the relationship between network topology and network controllability [10][11][12][13]. Researchers have proposed that all hub nodes with a high degree or betweenness centrality could be chosen as driver nodes [14]. Jalili et al. [15] found that an optimum driver node could not always be a hub node. To elucidate the configuration of driver nodes for the optimum network pinning control, a differential evolution method was used. The method worked well, but it was only suitable for undirected networks. Assuming that the objective network had finite-dimensional linear dynamics [16], the network could be structurally controlled by applying one time-varying input to the power dominating set. In the practical applications with economical and physical constraints, driver nodes can not always be freely selected to be injected to network nodes. In this context, Lo et al. [17] addressed a geometrical framework for the partial controllability issue of networks by solving an integer linear programme. The approach was also suitable to optimize the complete controllability of networks. The network permeability index provided a quantitative understanding of the challenge of controlling a network partially or completely.
The dynamics of the network could be expressed as a linear time-invariant system _ xðtÞ ¼ AxðtÞ þ BuðtÞ, where x(t) is the state vector of the network. Assuming a network has N state nodes and P control nodes, AR N×N is the coupled matrix between state nodes, u(t) is the control or input vector forced on the network, and BR N×P is the input matrix. The general approach for the controllability problem _ xðtÞ ¼ AxðtÞ þ BuðtÞ is to determine a proper input matrix based on the Kalman rank condition such that the pair (A, B) is controllable [10]. However, this controllability problem has a large computational load with 2 N possibilities assuming each node can be either driven or not driven [11], and this exponential growth is especially rapid when the network size is large. To overcome this difficulty, Liu et al. [10] introduced the structural controllability concept [18], which ensured that the Kalman rank condition was verified. They first found that the number of driver nodes for the full controllability of a complex network mainly depended on the network degree distribution. The process controllability of network dynamics was explored by transforming node dynamics into edge switch dynamics [19] and resulted in similar controllability conclusions to those obtained by Liu et al. [10].
The structural controllability methods based on graphical analysis of pair (A, B) for the system _ xðtÞ ¼ AxðtÞ þ BuðtÞ [18] could identify n D for arbitrary directed networks [10]. Several effective methods have been proposed to identify the minimum number of driver (control) nodes (n D ), for example, the maximum matching (MM) method [20], the cavity method [21], and an extremal optimization (EO) algorithm [22]. The computational load of determining n D could be effectively reduced based on the MM method, which has the computational complexity of Oð ffiffiffiffi N p LÞ, where L is the number of linked edges between state nodes. EO [22] was proposed based on the Kalman rank condition to identify n D for the full controllability of directed networks with the computational complexity of O(N 4 P 3 ). However, this structural controllability framework [10] is not applicable to undirected networks for the symmetric characteristic of the network matrix or networks with exact link weights [11,23,24]. These limitations prompt the development of exact network controllability theory, which is an exact controllability framework for the controllability of complex networks with arbitrary network structures and link weights. It optimizes the complete controllability of networks based on the Popov-Belevitch-Hautus (PBH) rank condition [25], which is an alternative criterion that is equivalent to the Kalman rank condition [26]. The PBH controllability method requires the sequential computation of the eigenvalues of the N × N matrix A and the rank of the N × (N + P) PBH matrix. The computational complexity of the eigenvalue computation of matrix A and the PBH matrix rank is O(N 3 ), O((N + P) 3 ), respectively [27]. Thus, the computational complexity of the PBH method is O((N + P) 3 ). The Kalman controllability method does not require an eigenvalue computation. However, it requires the rank computation of the N × NP Kalman matrix with the computational complexity of O(N 3 P 3 ), which is larger than that of the PBH controllability method. Representative exact structural controllability methods consist of a maximum multiplicity theory (MMT) [11] and an effective self-adaptive genetic algorithm (GA) [28]. n D was computed based on the MMT to be equal to the maximum geometric multiplicity of all eigenvalues of the network [11], and the computational complexity is O (N 2 (logN) 2 ). The GA [28] was studied to identify n D of arbitrary networks, and the computational complexity is O(2m × (N + P) 3 × l), where 2m is the population size and l is the number of different eigenvalues of the controlled network.
The authors demonstrated the evolution process of network topology of two networks (Erdős-Rényi (ER) and scale free (SF)) [28], i.e., the dynamic change of network topology with different number of control nodes being injected to network state nodes. The GA algorithm [28] exceeded the EO [22] much in terms of the convergence speed and iterations. However, networks with the more complexity than the scale 150 were not studied in [28] and the maximum network scale that EO processed was 200 [22]. And the convergence speed and iterations were still not satisfactory, for example, for ER, with both 150 state nodes and 150 control nodes, and the average degree 5.0, n D converged at the 101st generation after 398.03 s [28]. The results showed that it remained a challenge for the two algorithms to optimize networks with hundreds of thousands or even larger networks. Additionally, almost all real networks have small-world (SW) properties with a large cluster coefficient and short average distance [29] (e.g., power grids, transportation networks, and social networks). The addition of the SW network controllability study is also significant for better mimicking reality.
Therefore, based on the PBH rank condition, we propose a parallel quantum genetic algorithm (PAQGA) to more rapidly determine the minimum number of control nodes. The proposed algorithm is suitable for arbitrary networks that comprise both control nodes and state nodes. The simulation results for a series of benchmark networks demonstrate the effectiveness of the algorithm. Furthermore, we demonstrate the relationship between the controllability of a network and its network properties such as network average degree, degree heterogeneity, power-law index, and clustering coefficient.
The remainder of this paper is organized as follows: In Section 2, we provide a description of the issue in which a network can be controlled through a small amount of control nodes. In Section 3, we introduce the PAQGA for the solution of the minimum number of control nodes to exactly control arbitrary networks. In Section 4, we analyze and discuss the performance and experimental results of the proposed framework by studying popular ER, SW, SF, and some real-life networks. We draw conclusions and suggest future work in Section 5.

Problem definition
In this paper, we provide a descriptive definition of the entire controllability problem of a directed weighted network.
Definition 2.1 [22]. A network that contains P control nodes and N state nodes can be expressed as a triple tuple G = (V,E,W), where V = V s [ V c , V s = {v 1 ,v 2 ,. . .,v N } = {x 1 ,x 2 ,. . .,x N } is the set of state nodes, and V c = {v N+1 ,v N+2 ,. . .,v N+P } = {u 1 ,u 2 ,. . .,u P } is the set of control nodes; E = E s [ E c , E s 2 V s × V s is the set of the linked edges between state nodes, and E c 2 V c × V s is the set of linked edges between control nodes and state nodes, where each state node can only be connected to one control node; and W 2 R (N + P) × (N + P) is the set of edge weights, w ij = 0 if there is not a link between v i and v j ; otherwise, w ij (w ji ) represents the strength that Remark 2.1 [22]. The set W can be expressed by the representation of a block matrix that contains A and B as follows: Definition 2.2 [10,22]. The network G = (V,E,W) can be represented by where x(t) = (x 1 (t),x 2 (t),. . .,x N (t)) T is the state vector of the network, AR N×N is the coupled matrix between state nodes, u(t) = (u 1 (t),u 2 (t),. . .,u P (t)) T is the control or input vector forced on the network, BR N×P is the input matrix, and B = {b ij }, b ij is the weight of a directed link that the input signal u j (j = 1,2,. . .,P) points to the network state node x i (i = 1,2,. . .,N). For simplicity, hereafter the time symbol (t) will be omitted.  [22]. A control scheme D of a network G = (V,E,W) is determined by the selected control nodes with definite number and their acting position. D could be represented by a binary diagonal matrix as D = diag{d 1 ,d 2 ,. . .,d P }, where d i (i = 1,2,. . .,P) is a variable of value zero or one, and d i = 1 means that the control node u i is chosen to be a component of the network control strategy; otherwise, u i is removed together with its associated links.
Remark 2.2 [22]. Based on Definition 2.3, a novel control scheme D Ã is determined for which a different set of control nodes is selected. Then a novel network topology is generated where r is the number of selected control nodes. Accordingly, the network dynamics are also changed as where B Ã 2 R N×r is the new input matrix that represents the connections between new chosen control nodes and network state nodes, and u Ã 2 R r is a time-variable input vector that contains r control nodes. Definition 2.4 [22]. MR P×r is the index set of the selected control nodes, M = {m ij }, and m ij = 1 means that the j th chosen control node is u i , i = 1,2,. . .,P, j = 1,2,. . .,r, r P. Remark 2.3 [22]. M is constructed by the nonzero columns of the control scheme D Ã . For example, if u Ã = {u 1 ,u 2 ,u 3 } is chosen from a previous control node set u = {u 1 ,u 2 ,u 3 , u 4 } to be a new control scheme D Ã = diag{1,1,1,0}, then M is obtained from this D Ã as M ¼ ; ð4Þ Remark 2.4 [22]. The evolving input matrix and input vector can be revised as B Ã = BD Ã M and u Ã = M T u, respectively. Then Eq (3) can be rewritten as is satisfied for each different eigenvalue λ i of the state matrix A, where I N 2 R N×N is an identity matrix.
For an arbitrary network G = (V,E,W), our purpose is to determine the minimum control nodes to guarantee its full control. Based on the above analysis, the controllability problem can be transformed into a single target restricted optimization problem as subject to where A is the state matrix, B is the original input matrix, D is the original control scheme, M is the indicator matrix that is derived from the nonzero column of D, N and P are the dimensions of A and B, respectively, λ i is the eigenvalue that belongs to A, eig(A) is the set of different eigenvalues of A, d j is the element of D, and d j = 1 when u j is selected and d j = 0 otherwise.

Solution framework
Overview of the QGA GA is a global optimization method that can optimize problems with multiple parameters to reach near the global optima [30][31][32][33]. However, in some practical applications, it often requires multiple iterations because of the slow convergence speed and prematurity features, and easily falls into the local minima [34,35]. Additionally, for many complex problems, a large population is required to obtain the optimal solution. The convergence of GA mainly depends on the selecting operation, which largely affects the convergence speed. Additionally, its searching capability mainly relies on crossover and mutation operations, which primarily affect the occurrence of the premature phenomenon. Therefore, regarding enhancing GA search performance, the approach used to choose suitable selecting, crossover, and mutation strategies has been always an urgent and pivotal issue in the study and application of GA [33,36]. For small and medium-sized applications, the solution could be achieved within a tolerance range using GA. However, a gene (typically encoded with a 0-1 string) in a GA chromosome typically delivers certain information, which limits the population diversity. It performs worse in multivariate issues, for example, the controllability study of complex networks, which mostly has complex structures, and large-size nodes and links.
Combining quantum computing and GA, and adopting qubits as the representation of chromosome genes [37], QGA is a proper intelligent optimization algorithm for solving the network controllability problem [38]. These QGA qubits cover all possibilities for the linear superposition property of quantum information, which could reduce the algorithm's complexity and promote the achievement of the optimal solution under a smaller population [37].
In quantum computing, |0i and |1i signify two basic states of microscopic particles. According to the principle of the superposition property, the superposition state of quantum information could be the linear combination of the two basic states [39], which can be written as where α and β are the state probability amplitudes of a qubit, and α 2 and β 2 are the probability that a qubit changes to be state |0i and state |1i, respectively. One qubit also can be expressed Assume the number of optimization variables is n and the population size is 2m. The i th chromosome is denoted by G i (i = 1,2,. . .,m) as where . . . ; m. G i contains two parallel gene chains or individuals (α i1 ,α i2 ,. . .,α in and β i1 ,β i2 ,. . .,β in ). Each individual is a candidate solution of an optimization problem: ( QGA and enhanced QGA have already been studied to optimize many combinational problems [40][41][42]. For example, QGA overmatches classic GA with less complexity and higher performance in 0-1 combinational optimization problems [39]. Adaptive QGA models were proposed and tested on classical combinational problems, such as knapsack, maxcut and onemax [38], the multi-aircraft cooperative target allocation problem, and constrained engineering design problems [43]. However, the time efficiency was not seriously stressed. To increase the speed, a parallel QGA was developed and effectively applied to a knapsack problem [44]. It divided the entire population into subpopulations on different parallel processors and used the migration rate for the information exchange of these subpopulations. However, the Q-gate rotation was implemented according to a fixed lookup table, which did not take full advantage of the dynamic differences between individuals during the iterating process. Inspired by current achievements, to quickly and efficiently solve the controllability problem of complex networks, we investigated a PAQGA scheme, in which: 1) partial programs of the algorithm are executed in parallel; 2) a set of adaptive Q-gate rotation rules are proposed and adaptive crossover operation are used; and 3) population catastrophe is implemented to accelerate convergence.

Workflow of the PAQGA for network controllability
Based on the above, each control scheme D is a diagonal matrix, whose elements on the primary diagonal are either zero or one. Therefore, we adopt the binary mechanism to encode the algorithm chromosome G i (i = 1,2,. . .,m), and each binary gene chain can represent a control scheme, as shown in  To apply PAQGA conveniently, a penalty term Pen i (D) is defined to convert the optimization problem described in Eqs (8), (9) and (10) into an unconstrained optimization problem. Pen i (D) is used to evaluate the i th perturbation for a specific control scheme D: where σ i is the penalty coefficient defined as is a strictly increasing positive sequence to reduce the calculation burden of minimizing the penalty function, and l is the total number of distinct eigenvalues of A. The overall penalty of D is the sum of Pen i (D) expressed as According to the PBH rank condition, when the network G = (V,E,W) is fully controllable, Pen(D) should be zero and vice versa. Therefore, the fitness function can be achieved by merging the penalty term into the optimization Eqs (8), (9) and (10): For the optimization problem, our objective is to minimize f(D) based on the 0-1 integer values of d j (j = 1,2,. . .,P). After defining the chromosome representation and fitness function, the network controllability problem can be optimized using the following steps. Fig 6 shows the flow chart of the proposed PAQGA for the optimization problem.
Step 1: The population at the t th generation is denoted as QðtÞ ¼ fG t 1 ; G t 2 ; . . . ; G t m g; where N is the number of qubits, that is, the number of network state nodes, and maxgen is the maximum iterating generation. Initialize the initial population as All the quantum states (α ik and β ik ) in the PAQGA are initialized in parallel with the value 1 ffi ffi 2 p , i = 1,2,. . .,m, k = 1,2,. . .,N. Additionally, set D best = D 0 , the iterative generation t = 1, and σ 1 = 10P.
When we set the initial control scheme D 0 = diag{d j = 1},j = 1,2,. . .,P, the initial best fitness value is f(D best ) = P + 0 = P. It is easily proved that with σ i = {cσ i−1 } (σ 1 = 10P, c > 0), the fitness f(D) P if and only if the control scheme D always satisfies the constraint Eq (9). With the initialization f(D best ) = P, whenever D best is updated by D i , we have f(D i ) < f(D best ) = P, which means D i meets the constraint Eq (9). Thus, D best always evolves in the feasible region that makes the network entirely controllable.
Step 2: Observe the qubits of each individual of Q(t 0 ) in parallel following the rules in Section 3.3 and obtain the binary strings.
Step 3: Evaluate each individual of Q(t 0 ) in parallel and save the optimal individual as the evolving goal in the next generation.
Step 4: While (maxgen is not reached), do the following: (a) Observe the qubit value of each individual of Q(t) in parallel following the rules in Section 3.3.
(b) To increase the diversity of the population and inherit the excellent genes from the previous population, the adaptive crossover operation is performed in parallel in accordance with Section 3.4.
(c) Evaluate each individual of Q(t) in parallel and store the optimal individual as the evolving goal in the next generation.
(d) Perform the Q-gate rotation operation in parallel and obtain the offspring population Q (t+1).
For each individual, two parallel gene chains update simultaneously. Rotation angle θ is first computed based on Section 3.5 and then the qubits are updated by Q-gate rotation. The Q-gate is expressed as [45] QÀ gate ¼ cosy À siny siny cosy The Q-gate rotation operation is

Observing operation
Each qubit of the chromosome can be adjusted to be at a stationary state using an observation operation. We adopt a random observing method by running the following pseudocode in parallel: where rand(k) is a random digit. If rand(k) is not less than ðb t ik Þ 2 (the probability to be state |1i), then the observed value of the qubit a t ik is 1 and 0 otherwise.

Crossover operation
The crossover operator is an important operation of GA. Information about individuals can be exchanged using the operation. Subsequently, excellent genes could be reserved for population evolution to move in a better direction. To increase the diversity of the population and improve the optimization performance of PAQGA, the crossover operator is introduced. We obtain novel binary values a 0 ik by crossing each binary qubit value α ik ,i = 1,2,. . .,m, k = 1,2,. . .N with corresponding information on the historically best control scheme, that is, D best (k,k), based on a certain crossover probability. The specific crossover mode is In the early days of population evolution, there existed relatively big differences between individuals. Therefore, the crossover possibility to produce better offspring should have been bigger. Moreover, if we increased the crossover probability at this time, the evolution process would have been accelerated. By contrast, in the late stages of evolution, differences between individuals became smaller as the best solution was approaching. The crossover probability should have been correspondingly diminished to reserve the good genes. We design an adaptive crossover operator as where p c (i) is the crossover probability of the i th current individual, Q i is the number of those individuals whose fitness is better than that of the historically best individual, p c0 is the initial crossover probability, f max and f min are the previous worst fitness and best fitness, respectively, and f(G i ) is the fitness of the i th current individual. We can observe that p c (i) becomes bigger when the control scheme G i becomes worse and vice versa. Moreover, p c (i) is inversely proportional to Q i , which means that if there are not so many good individuals, p c (i) should be bigger to produce a greater number of better individuals; otherwise, it should be smaller because the evolving individuals are becoming better. The improved adaptive crossover operation from Eq (20)  ( A better population is determined after the crossover operation following the pseudocode.

Rotation angle updating rules
Learning from the solid lookup rules [39], we present a set of adaptive rotation angle updating rules in Table 1. The rotation angle θ i (θ i = δ i Ã s(α i ,β i )), i = 1,2,. . .,m dynamically varies according to the evolution process.
If D c (i) 6 ¼ D best (i), the rotation angle δ i is adaptively proportional to f ðD c Þ f ðD best Þ . If f(D c ) < f (D best ), the angle will be smaller; otherwise, it will be bigger. Initially, a big initial angle is set.   As the iteration proceeds, the differences between individuals decrease and δ i becomes smaller. In this way, the probability amplitude evolves in the direction of the optimal solution.

Population catastrophe
When the best individuals in several successive generations are identical, it shows that the algorithm falls into a local minimum. At this moment, catastrophe operations for the current population should be performed to take it out of the constraint and start a new search. Specifically, the successive best individual is retained in the new population Q(t + 1) and the remaining individuals in Q(t + 1) are regenerated as a large disturbance. The pseudocode of the catastrophe operation is as follows: start obtain the best individual corresponding to the optimal fitness; keep this best individual; rebuild the rest in parallel; end The strategy would prefer that the population eliminate its dull state rather than make it degenerate, which is an effective means to commence a new search.

Simulations and analyses
We used the orthodox ER random [46], SF [47], SW networks of NW type [48] and some realworld networks as benchmarks to illustrate the feasibility of the PAQGA for optimizing the controllability of arbitrary networks that encompass control nodes and state nodes. Additionally, we also conducted an analysis of the relationship between the network topology and number of control nodes. ER and SF networks were obtained from the static model [49] with N state nodes and P candidate control nodes (N = P). Each control node pointed to state nodes with uniform probability and the weights of all edges were randomized between zero and one. SW networks were generated from randomized adding edges [48,49]. The characteristics of random regular networks, ER, SF, and SW networks are illustrated in Fig 8. We define the number of selected control nodes that correspond to the current best control scheme as n c and the density of these selected control nodes as N c , where N c = n c /N. The minimum number of selected control nodes after the optimization process is denoted as n cm , and the minimum control node density is N cm = n cm /N. To implement the parallel strategy, we performed the following simulations on eight MATLAB 1 workers. The parameters of the PAQGA were set to 2m = 30, maxgen = 100, and p c0 = 0.06. All the following experimental results are the average of 10 independent simulations and the standard deviation is 0.01.

Performance of the PAQGA
To show that the optimal solution (D best ) at each generation always evolves in the feasible region, that is, Pen i (D best ) = 0, we conducted experiments on different networks. All these networks were directed with 100 state nodes and 100 candidate control nodes. The experimental results are shown in Fig 9. The figure shows that the best fitness quickly converged to the minimum value after approximately the first few generations. The mean current fitness fluctuated dramatically because of the operations of qubit cross, qubit catastrophe, and Q-gate rotation. The penalty was always equal to zero, which means that D best always met the PBH rank condition throughout the entire optimization process.
We compare the performance of PAQGA of optimizing network controllability with that of EO [22] and adaptive GA [28] on a list of popular networks and real-life networks. The comparison results are tabulated in Table 2.
From the columns of n cm , it can be observed that the three algorithms almost converged to the same value, which demonstrates that PAQGA, GA, and EO all had a good ability to determine the optimal control nodes. When the size of the network was small (e.g., N = P 50),  PAQGA took slightly more time than GA and EO to determine the best solution. This paradoxical phenomenon is attributed to the launching of the MATLAB 1 distributed server, and the launching time was approximately 2s. However, once the server started, PAQGA showed a greater advantage in processing large-size networks over GA and EO. For example, for the ER network with hki = 6.0 and N = P = 200, PAQGA obtained n cm at the fifth generation and took 12.49s; for the same network, GA took 873.05s at the 128 th generation and, and EO required 29 generations and 1545.84 s. By comparing the computational time, PAQGA saved 98.57% more than GA and 99.19% more than EO.
A parallel version of GA was transformed from the adaptive GA [28] using η MATLAB1 workers with the computational complexity of O 2mÂðNþPÞ 3 Âl Z , where 2m is the population size, l is the number of different eigenvalues of the controlled network. We compared it with the proposed PAQGA, and the results are shown in Table 3. It can be inferred from Table 3 that the parallel computation (allowing for multiple processors) contributes to the performance of algorithms. However, it is not the only important factor. The computational efficiency of EO, GA, parallel GA and PAQGA could be reflected by their computation complexity. First, the computation of the PBH rank matrix in GA, parallel GA and PAGQA and the Kalman rank matrix in EO is the most time-consuming. This is the main factor affecting their speedability. The rank computation of Kalman matrix takes much Table 2. Performance comparison of PAQGA, GA, and EO on different networks in terms of n cm , the minimum iterating generations, and computational time.
Power-law index of SF networks in these experiments was γ = 2.1. '/' indicates that the corresponding results were not available for the computational time limit. For data sources, see Supplementary information S1 Dataset.  Table 4.
It can be seen that N cm of PAQGA agrees with that of MMT on these real-world networks, and is slightly greater than or equal to that of MM. The experimental results show the efficiency of PAQGA in identifying the minimum control nodes. Nevertheless, PAQGA is at a disadvantage in computational time compared with MMT and MM, although such defect could be improved by adding the number of processors or using computer groups. For example, the cost of calculating N cm of Wiki-vote network is reduced to 301.34s with 12 processors. PAQGA is an intelligent probabilistic optimization algorithm that provides approximate solutions. The optimal solution cannot be guaranteed to be found. Moreover, the optimal solution of a problem is typically unknown in advance. We used the structural controllability theory [10,18] as the benchmark to compute n cm to test the validation of PAQGA. The results are shown in Fig 10. We can observe that for ER, SF, and SW, the obtained n cm is the same as the benchmark result, which indicates that the proposed PAQGA was effective in determining the minimum control nodes of complex networks.

Discussion and analysis of results
Applying PAQGA, the optimization results and evolution process of network topology can be achieved. The results are intuitively displayed in Fig 11. From Fig 11(A), N c of different networks quickly converged to a steady minimum value, which indicates the effectiveness of the PAQGA. For example, N c of ER with <k> = 4.0 converged to 0.05 at the fifth generation, which demonstrates that five control nodes were sufficient to maintain network controllability. For SF with <k> = 6.0 and γ = 2.1, N c rapidly decreased to a minimum value of 0.07 at the seventh generation, which means that at least seven control nodes were required to fully control the network. Fig 11(B)-11(D) together capture the evolution of the SW network with <k> = 6.0 at the zeroth, fifth, and seventh iteration, and the convergence trend of the control nodes can be acquired.
We also found that two networks with different <k> required different N cm . For example, for ER with <k> = 4.0 and <k> = 6.0, N cm of the network with <k> = 4.0 was 0.05 and with <k> = 6.0, N cm = 0.02. Additionally, N cm of SW networks with <k> = 4.0 and <k> = 6.0 was 0.25 and 0.22, respectively. Second, for networks with the same <k> and different γ, N cm also differed. For example, consider SF with same <k> = 6.0, and different γ = 2.1 and γ = 3.0. The two networks had N cm = 0.06 and N cm = 0.04, respectively. Third, for networks with the same γ and different <k>, N cm was also different, which can be determined from Table 2. These results led us to conjecture that N cm had a relationship with <k> and γ.
To confirm our hypothesis, we performed simulations on a set of different networks and plotted N cm as the function of <k> and γ. The results are shown in Fig 12. From Fig 12(A), it is obvious that N cm of networks with fixed γ decreased monotonically with <k> until N cm became slowly flat. Additionally, the downward trend was of asymptotic exponential dependence, which suggests that the sparse network required more control nodes to maintain full controllability. From Fig 12(B), we can observe that N cm with fixed <k> decreased as γ increased. The results indicate that N cm may be influenced by the degree heterogeneity, denoted by H, which is the standard deviation of the network node degree distribution [57]. In this paper, H is defined as where k i is the degree of state node i.
To determine the relationship between N cm and H, we examined N cm as a function of H and obtained the results shown in Fig 13. From Fig 13(A), it can be observed that a larger N cm always corresponded to a larger H and smaller γ. Fig 13(B) shows that the network with a smaller hki and larger H typically required a larger N cm . The results suggest that the larger the differences between node degrees, the more control nodes were required to entirely control the network.
SW networks have the remarkable characteristics of a large clustering coefficient, denoted by C, which represents the overlapping degree of friend circles of two adjacent state nodes and is defined as where i is node i, k i is the number of edges between node i and other nodes, and E i is the number of edges among the k i nodes. For SW networks, N cm may be affected by C. To explore the relationship between N cm and C, we plot N cm as a function of hki and C,shown in Fig 14. From Fig 14(A) and 14(B), we can determine that a larger C corresponded to a smaller N cm , which indicates that the more interconnected the network, the fewer control nodes were required to control the network. For other networks, such as ER, SF, the conclusion also holds.
Considering the aforementioned analysis results together, we can determine that for a given network with both control nodes and state nodes: 1) the sparser the network, the more control nodes were required to control it; and 2) the more heterogeneous the network, the more control nodes were required to guarantee its full control. We reflect that the sparse and heterogeneous network is the most difficult for guiding its dynamic evolution (see Tables 2  and 4 and Figs 10(A) and 12(B)). The consistency between the results from our approach and from these existing methods [10,11,22,28] confirms the similarity between them for directed networks, which not only further validate these existing methods, but also reflect the effectiveness of our method.
To evaluate the controllability of directed networks, the structural controllability framework based on the MM method is still the best for its error-free feature [11]. Like the MMT, the PAQGA also relies on the eigenvalues and the rank of the network matrix, the computation of which inevitably introduces numerical errors. Further, MM and MMT both surpass PAQGA in computational efficiency in identifying both the minimum set of driver nodes and the number of these driver nodes. However, the PAQGA can have a wider range of applications. For example, the PAQGA is valid for networks containing a number of self-loops with identical or different weights, and networks with bidirectional connections between two nodes. The PAQGA is also applicable to undirected networks, where the structural matrix assumption is slightly violated because of the network symmetry. Further, combined with advantages of computer hardware and the adaptive strategies itself, PAQGA has great room for improvement. Taken together, the PAQGA as an alternative exact structural controllability framework provides us deeper understanding of the controllability of complex networked systems.

Conclusions
In this paper, we introduced a PAQGA to optimize the controllability of arbitrary networks with control nodes and state nodes under the PBH rank condition. In addition to MATLAB1 workers, more parallel mechanisms can be flexibly embedded in the PAQGA, for which more threads concurrently processing could further promote the time efficiency of generating a solution. Analyses and simulation comparisons demonstrated the effectiveness and applicability of the proposed PAQGA. Furthermore, we found that the minimum control nodes were affected by the network degree distribution, degree heterogeneity, and clustering coefficient. The sparse and heterogeneous network is the most difficult to be fully controlled.
In our study, the topology that comprises state nodes remained static during the entire evolution process. However, networks normally evolve over time, which manifests as the increasing or decreasing of different nodes and their links. In the future, we will focus on the controllability of dynamic networks. Furthermore, we hope to explore how to use the obtained minimum control nodes to steer an intermediate network to evolve into our desired network considering realistic energy constraints.