Optimizing target nodes selection for the control energy of directed complex networks

The energy needed in controlling a complex network is a problem of practical importance. Recent works have focused on the reduction of control energy either via strategic placement of driver nodes, or by decreasing the cardinality of nodes to be controlled. However, optimizing control energy with respect to target nodes selection has yet been considered. In this work, we propose an iterative method based on Stiefel manifold optimization of selectable target node matrix to reduce control energy. We derive the matrix derivative gradient needed for the search algorithm in a general way, and search for target nodes which result in reduced control energy, assuming that driver nodes placement is fixed. Our findings reveal that the control energy is optimal when the path distances from driver nodes to target nodes are minimized. We corroborate our algorithm with extensive simulations on elementary network topologies, random and scale-free networks, as well as various real networks. The simulation results show that the control energy found using our algorithm outperforms heuristic selection strategies for choosing target nodes by a few orders of magnitude. Our work may be applicable to opinion networks, where one is interested in identifying the optimal group of individuals that the driver nodes can influence.

of control signal increases 18 . Thus, attaching additional control signals onto a networked system beyond the minimum driver node set is one way to achieve network control with reduced energy cost.
The way in which the additional control signals are attached can also result in reduced control energy. Lindmark and Altafini considered the eigenvalues of the network and proposed strategies for selecting the placement of additional control signals to minimize control energy cost 20 . Chen et al. analyzed stems, obtained from minimum driver node set 6 , and calculated all possible direct shortest paths from driver nodes to non-driver nodes to obtain the Longest Control Chain (LCC), which is the longest path from the obtained all possible shortest direct paths 21 . They found that by adding additional control signals in such a way that the length of the LCC is minimized, control energy can be reduced significantly. Li et al. proposed an algorithm using matrix derivative projected gradient descent to iteratively search for the energy optimal placement of control signals 17 . This optimization model was later simplified by Ding et al. 15 and Li et al. 16 . In 2018, Li et al. proposed an improved and generalized approach based on previous works to again obtain the energy-optimal placement of control signals 22 . The problem of reducing control energy by strategic placement of additional control signals has been extensively researched. However, in all of these works, complete control was considered, where the control signals steer the full node set towards the predefined goal state vector.
While full control may be necessary in some types of engineered systems 23 , controlling just a subset of nodes (typically termed target control or targeted control) may be more sensible in large complex dynamical systems. In 2014, Gao et al. proposed an alternate k-walk theory and a greedy algorithm to obtain the minimum number of driver nodes to control just a subset of the full node set in a complex network 24 . In 2015, Iudice et al. presented a geometric framework to find the driver node set, limited by practical constraints, to reach as many target nodes as possible 25 . In 2017, Liu et al. considered the controllability of the giant connected component in a directed complex network 26 . With regards to energetic considerations, Klickstein et al. showed that the energy cost scales exponentially with the cardinality of the target node set 19 . Gao et al. proposed an algorithm to obtain the placement of control signals to optimize the energy cost when controlling just a subset of nodes in directed complex networks 27 . Furthermore, recently, the controllability Gramian of lattice graphs was studied, which turns out to be useful when the number of target nodes is small 28 . While target control has been researched extensively, how to optimize control energy with respect to the selection of target nodes in a complex network is still an open question.
Thus, different from all previous works, we wish to ask a slightly different question: How can we pick the target node set in such a way that the control energy is minimized? To that end, we will be employing a cost function optimization model based on projected gradient descent, similar to earlier works [15][16][17]22,27 . The main difference between our work and the existing literature is the variable matrix being optimized: While previous works have focused on optimizing the choice of input signal (matrix B), we optimize the choice of target nodes (matrix C). Furthermore, previous derivations of the index notation gradient information uses the I-Chain rule, which can be difficult to understand. Here, we derive the index notation gradient information using the standard chain rule and product rule in a general way, which is simpler to follow.
In this paper, we will be examining how the choice of target nodes could be optimized to minimize the energy cost function. Using the formulated energy cost function, we derive the matrix derivative of the energy cost function with respect to matrix C, which is the energy cost function gradient information. With the gradient information obtained, we perform an iterative search using the trace-constraint-based projected gradient method (TPGM) proposed in Ref. 22 to obtain target nodes which are energetically favorable. We then compare the energy cost that we would have gotten from choosing to target control nodes using a heuristic selection scheme such as random selection, and node degree-based selection. Our simulation results show that the solution obtained from TPGM reduces the energy cost by a few orders of magnitude.

Problem formulation
In standard complex network controllability literature, we are interested in studying N coupled system of equations whose dynamics are linear and time invariant (LTI): where x(t)∈R N×1 is the time-varying state vector, y(t)∈R P×1 ( P ≤ N ) is the the subset of time-varying output state vector that we want to target control. u(t)∈R M×1 ( M ≤ N ) is a time-varying external control signal which we use to drive the network. The time invariant A ∈ R N×N matrix represents the network topology, where for directed complex network we have nonzero matrix element {a ij } when there is a directed link from node j to node i, and zero element if no link exists from node j to node i. Since many complex systems tend to enjoy passive stability, our modeling also allows for self-links, where {a ii } is a negative real number 12 . The time invariant matrix B ∈ R N×M reflects the coupling between nodes and M number of external control signals u(t) , where matrix element {b ij } = 1 if node i receives a time-varying control signal from control input j, and zero element otherwise. Finally, we have the the control matrix C ∈ R P×N which tells us the choice of the P number of subset of nodes that we wish to target control, where matrix element {c ij } = 1 if node j is the i-th node out of all possible P nodes that we want to target control and zero otherwise. We note that matrix B and matrix C are column and row linear independent, where if {b ij } = 1 ( {c ij } = 1 ), then no other nonzero entries may exist for column j and row i in matrix B (C). We require this linear independence to reflect our modeling choice for one-to-one connection between control signals and nodes, and for one-to-one connection to individually target control each of the P number of nodes 19 .
We are interested in the energy cost needed to drive the states of the network using the control signal u(t), which is defined to be the cost function 9 (1)

Scientific Reports
| (2020) 10:18112 | https://doi.org/10.1038/s41598-020-75101-w www.nature.com/scientificreports/ This is the energy cost function that we want to optimize, with respect to the choice of target nodes (matrix C). Similar to previous works which considered trace constraint projected matrix gradient cost function optimization [15][16][17]22,27 , here we assume that the initial state vector x 0 ∼ N(0, 1) . The energy cost is where t f is the total control time and initial time is set to zero throughout the rest of the paper, t 0 = 0 . The energy cost is a function of the control trajectory, i.e. the state space pathway between x 0 and y f . It should be noted that (3) is the expected energy cost over all realizations of initial state vector x 0 , picked under standard normal distribution. In general, the cost function E (t f , B, C) is dependent on the matrix B and the final time t f . However, in this paper, we are focusing on optimizing the cost function with respect to matrix C while keeping input matrix B and final time t f fixed. Thus, for clarity, we drop these two dependencies and write the cost function only as a function of either E (C) or E (C T ) throughout the paper. Consequently, the energy-optimal control signal, derived from optimal control theory 10,19 can be calculated as The N × N controllability Gramian matrix is an important quantity in control theory, well-known to be real, symmetric, and semi-positive definite. In practice, we can use the controllability Gramian matrix to verify that the choice of matrix B is ensuring network (target) controllability by checking that the (output (CWC T ) ) controllability Gramian matrix is invertible 9,19 . We optimize the cost function as follows: where we have constrained the solution produced by the gradient descent iterative algorithm to stay on the manifold tr(CC T ) = P.

Results
The energy cost function can be obtained by substituting the energy-optimal control signal u * (t) into the cost By varying Eq. (7) with respect to C T , we find that The full derivation can be found in the Supplementary Information. Using (7) and (8) and applying it to the TPGM algorithm, we solve for the energy-optimal target node matrix and obtain the optimal solution, C * . In the simulations, we have set the control time to be t f = 2 , and the desired final output state vector is Algorithm. To optimize the control energy by varying the selection of target nodes, we modified the traceconstraint-based projected gradient method (TPGM) formulated by Li et al. 22 , to focus on C T , the target nodes choice. Let C T be the basis of the target control matrix C T . We can obtain C T by performing Gram Schmidt orthogonalization on matrix C T29 . Define the projection operator TCT = (I N −C TC ) , where I N is the N × N identity matrix. The operator TCT projects any arbitrary matrix onto the space that is perpendicular to the manifold tr(CC T ) = P . For two arbitrary matrices A and B which have the same dimension, an angle between them can be defined as: where 0 ≤ θ ≤ π , . F denotes Frobenius norm, and we note that the matrices are perpendicular if θ = π/2 = 90 • . The TPGM optimization is given in Algorithm 1. www.nature.com/scientificreports/ In TPGM step 1, by random initialization of matrix C T 0 , we mean that we start with a N × P matrix of zeros and randomly set matrix element [C T ] ij = 1 , where node i is chosen to be the j-th target node. Once node i is chosen, we maintain row and column linear independence, and thus do not allow for the selection of the same node i to be target controlled. In the numerical experiments that we tried, while starting from a pure random dense matrix is possible when system size is small or the network topology is sparse, it is computationally inefficient. Furthermore, when starting from a pure random dense matrix, we are not able to obtain a sensible sparse optimal matrix at the end, when converting from dense optimal solution C * to sparse binary optimal solution (see "Selecting binary optimal target node set from C * " section).
In TPGM step 4, the numerical value of the learning rate or step size η is chosen empirically. When η is too large, convergence is not guaranteed and the algorithm will fail. While convergence is guaranteed when η is sufficiently small, the time taken to complete the iterative search will suffer if η is too small. Typically, we choose the learning rate to be between 1e−8 and 1e−3 , depending on the fraction P/N to be target controlled as well as the complex network topology. In general, the learning rate η can be varied to speed up the iterative process. For example, starting at η = 1e−8 , and then changing to η = 1e−4 when enough iterations have been run. As observed during experimentation, the convergence of η scales inversely proportional to the numerator of the . In TPGM step 4, Ĉ T k+1 is a non-normalized quantity, while in step 5, the updated solution is constrained onto the manifold surface tr(CC T ) = P , obtaining the normalized C T k+1 . In TPGM step 6, the angle θ k between k-th step gradient matrix ∇E (C T k ) , and the projected gradient matrix, TCT k ∇E (C T k ) is calculated to check for convergence. If cos(θ k ) approaches zero, or equivalently, when θ k approaches π/2 = 90 • , the algorithm is deemed to have converged. Numerically, the while loop terminating condition, ξ , refers to a small positive quantity, for example, ξ = 1e−2 . For some networks where the algorithm is unable to converge towards ξ = 1e−2 , we may relax this condition to ξ = 1e−1.
TPGM will iteratively update the initial proposed solution C 0 in the direction of quickest decreasing energy cost based on cost function derivative with respect to matrix variable C. When the search has finally converged, the obtained optimal target control matrix C * corresponds to reduced control cost of a dense matrix of real numbers, which corresponds to many-to-many connections from output nodes to complex network nodes of varying link strength. Based on the obtained optimal solution C * , the challenge is to obtain a one-to-one connection of output nodes and nodes with unity link strength, a sparse binary optimal solution, C * binary , while maintaining the characteristics of reduced control cost.
Selecting binary optimal target node set from C * . We propose two methods to find C * binary from C * . The first method is based on suppressing insignificant matrix elements in C * , and then evaluating the normalised quantity importance index vector, r: is the non-normalised importance score of node i being selected as a target node, and we are taking the row summation of the absolute of the transpose matrix, C T . |[C T ] ij | represents the numerical contribution towards C * for node i, target node j. Therefore, to find binary optimal solution that remains similar to C * , we want to find the nodes with the highest numerical contributions.
In our experimentation, we found that directly passing C * into (10) to find C * binary usually does not allow us to find the optimal target node set. This is because of the row summation of the absolute of the numerical contribution. For example, if row a has predominantly absolute values of around 0.2 in all its columns, while row b has an absolute value of around 0.9 in one of its columns, and mostly negligible numerical values close to zero in all of its other columns, then by Eq. (10), node a will be ranked higher than node b. Based on experimentation with elementary topologies, we find that numerical contribution characteristics similar to row b usually correspond to an optimal target node. Thus, we want the row summation to reflect that: we suppress the numerical values in each matrix element of C * if they fall below a certain value. (For a more detailed discussion, see Supplementary Information.).

Scientific Reports
| (2020) 10:18112 | https://doi.org/10.1038/s41598-020-75101-w www.nature.com/scientificreports/ We compute the suppressed optimal matrix by setting where σ is computed from the standard deviation of the absolute of all matrix elements [C * ] ij , and d = {0.0, 0.1, 0.2, . . . , 1.1, 1.2, . . . , 3.0} . For each suppressed matrix, we check that the rank of the matrix is P, before considering it as a viable candidate solution. If the suppressed matrix has rank P, then we pass the suppressed matrix into (10), and pick P nodes with the highest importance index. The case when k = 0 is similar to the proposed importance index vector formulation in Refs. 16,17,30 for finding optimal driver nodes.
The second method is based on selecting the matrix elements with the largest absolute numerical values. The process is as follows: First, begin with optimal matrix transpose, |[C * ] T | . Then, search for the absolute largest matrix element in each column, and order the columns in descending order, keeping track of the associate row indices. Start with a N × P zero matrix, C T binary , and set the matrix element for position (i, j) to be one, starting from the columns with the largest absolute matrix elements |[C * ] T ij | . At each step, check that for assigning position (i, j) to be one, no other nonzero element exists along row i or along column j. If there exists another nonzero matrix element along the row i or column j, then do not assign position (i, j) to be one, and record down column j which did not get filled. If rank (C T binary ) = P , then stop the process; otherwise, repeat the process described for unfilled columns j's for the case of second largest absolute matrix elements, third largest,...., until rank (C binary ) = P.
At the end, we pass all the obtained binary target node set candidate solutions into the objective cost function, Eq. (7), and pick the target node set which yields the lowest energy cost, denoting it [C * binary ] t , where t = {1, 2, . . . , 10} represents each independent iterative search. For the simulation results showing C * binary , we choose the best solution out of all [C * binary ] t .

Numerical experiments on elementary network topologies.
To understand the arrangement of energy optimal target node set, we perform the TPGM iterative search on elementary network topologies 6,30 .
They are: directed stem, circle, and dilation. A stem requires just a control signal, placed at the root, to become mathematically (but not necessarily numerically) controllable 11,26,31 . A circle is controllable with just one control signal attached to any of the nodes in the circle. A dilation requires a minimum of two driver nodes to become controllable. We select the driver nodes to be placed in such a way that their path distances are evenly spaced, which corresponds to the most energy efficient set up as the longest path distance from driver to non-driver nodes is minimized 21 . The set of driver nodes is assumed to be fixed (B fixed), and the target node set, consisting of 66.7% of all nodes, is our variable, e.g. find C.
For each of the 8 elementary topology labelled (a1-c2) to as shown in Fig. 1, we repeat the numerical experiment independently for 1000 times in order to find the optimal target node set C * binary , given that driver nodes position are fixed. (a1) is a N = 9 directed stem network, with M = 3 driver nodes, and P = 6 optimal target nodes. Evidently, as seen in Fig. 1, the energy optimal target nodes tend to minimize their geodesic 32 path distances from the driver nodes. Our findings are consistent with literature, which states that control energy cost increases exponentially with path distance 21,33 .
To test the robustness of the algorithm, we repeat the experiments for elementary topologies with fixed driver nodes placement similar to configurations (a1,a2,b1,b2,c1,c2) shown in Fig. 1, but with total target control nodes changed from P = 6 to P = 5 . In general, the conclusion that energy-optimal configurations correspond to minimized path distances from driver to target nodes still holds. For the stem network with 3 driver nodes, the optimal configurations of target nodes are {1, 2, 4, 5, 7} , {1, 2, 4, 7, 8} , and {1, 4, 5, 7, 8} , given in ascending order of energy costs, with the first node set being the true optimal. TPGM was able find the optimal configurations with probabilities 30.8% , 15.2% , and 0% respectively. For stem network with 2 driver nodes, the optimal configurations are {1, 2, 3, 5, 6} and {1, 2, 5, 6, 7} , with the former being the true optimal. TPGM success rates are 15.1% and 55.3% respectively. Thus, while control energy is minimized when target nodes are close to their nearest driver nodes, for stem network, there is a slight reduction in control energy when the target nodes are also closer to node 1, which is the root node needed to ensure controllability. It should be noted that the energy cost difference between local optimal configurations and global optimal configurations is typically smaller than 10% . For circle topology with 3 driver nodes, the optimal configurations are {1, 2, 5, 7, 8} , {2, 4, 5, 7, 8} , {1, 2, 4, 5, 8} and {1, 2, 4, 7, 8} , {1, 2, 4, 5, 7} , {1, 4, 5, 7, 8} , with the first three node sets being the true optimal. TPGM was able to find node sets belonging to true optimal energy cost 29.6% of the time, and optimal configurations with probability 5.7% . Circle topology with 2 driver nodes has optimal configurations {1, 2, 3, 5, 6} and {1, 2, 5, 6, 7} , with the former being the true optimal. Correspondingly, the search success rates were 5.5% and 34.0% respectively. While a circle is controllable with just one driver node placed anywhere, the difference in energy cost can be accounted for by grouping all output nodes to their nearest driver nodes: {1, 2, 3, 4} to node 1, and {5, 6, 7, 8, 9} to node 5. Thus, for energy-optimal control of 5 target nodes, the first 4 have to be picked to be close to their nearest driver nodes, and for the fifth node, a slight preference is given to picking target nodes belonging to the group with the lower directed path distance. In all of these, the control energy difference between true optimal and optimal are at most 2% . For dilation networks with 2 driver nodes, the optimal configurations are {1, 2, 6, 7, 8} and {1, 2, 3, 6, 7} , with the former being the true optimal. They were respectively found by the algorithm 32.   Simulations on a small random network. Next, we repeat the numerical experiment 1000 times for target controlling 5 nodes in a N = 10 random network (ER network), with average degree �k� = 2 , as shown in Fig. 2. The driver nodes position are fixed, and generically ensures full controllability of the network, regardless of the choice of C. The global energy-optimal configuration of target nodes placement is {1, 2, 4, 8, 10} and TPGM was able to find this configuration 7.1% of the time. However, TPGM was also able to find the next-best performing configuration, which yields energy cost of the same order of magnitude, {2, 4, 6, 8, 5} , with success rate of 13.9% . Thus, generally, it can be seen that when path distances are reduced, control energy is minimized. We quantify the performance of TPGM on the small random network in Fig. 3. In Fig. 3a, we display the probability mass function of TPGM obtaining all the solutions found, C * binary , and their associated energy costs, indexed in ascending order in energy, TPGM E index i . It should be noted that all 25 TPGM energy cost indexes are found from 1000 independent iterative searches, and are not exhaustive. Figure 3b shows the energy costs associated to each TPGM E index i . The full range of all possible energy costs is 10 C 5 = 252 , so 25 is roughly the top 10% best solutions. For comparison, the average energy cost of selecting 5 target nodes randomly without repetition, E(C rand ) , is also plotted and represented by the black dashed lined on the same graph. In Fig. 3c, we show all 252 energy costs found from a brute force search. Matching the solutions found from TPGM to the true list of all possible energy costs, we find that the top 6 solutions from TPGM belong to the true top 10 energy costs, as shown in Fig. 3d. Therefore, while the rate of TPGM finding suboptimal solutions is high, the energy costs performances of the suboptimal solutions are in general comparable to the true optimal energy cost, E index 1 . From inspection, we see that the top 15 solutions have very similar energy, which may explain why TPGM only found the global solution (index 1) 7.1% of the time. In addition, the solutions found from the search algorithm also tend to lie within the neighborhood of the lower end of all possible energy costs. Besides, they also outperform the random selection scheme by at least a few orders of magnitude.
Simulation results on complex networks. We apply the TPGM to different complex networks such as random networks (ER network), scale-free networks (SF network), as well as various real networks spanning a diverse range: electronic circuit networks, food web networks, and social networks. For each network, we perform the iterative search 10 times and compare the performance of the obtained solution C * as well as the associate binary optimal solution C * binary and compare them to heuristic selection schemes such as random selection (repeated over 100 independent realizations) and degree-based selection of target nodes. We select 40% of the nodes to be driver nodes using the standard way 6,7 to ensure controllability, and pick the remaining nodes randomly. For each network, once the driver nodes are picked, they remain fixed.
For the model networks of N=100 SF ( γ = 2.8 ) and ER, with average degrees k av = 2.5 , we plot the control energy needed for controlling target node set of varying cardinality P/N = {0.1, 0.2, . . . , 0.9, 1.0} in Fig. 4a,b. Consistent with the previous findings 19 , we find that the control energy scales exponentially with the cardinality of the target node set, regardless of the target node selection scheme chosen. However, when comparing the control energy of various selection strategies, our results show that generally, as compared to heuristic selection schemes, the energy cost for controlling target node set as found by TPGM, C * and C * binary , is lower by a few orders of magnitude. For the result of Fig. 4b, the performance of the in-degree descending selection scheme (meaning that we choose P/N% of nodes that are ranked within the top P/N% of largest in-degrees) in the 4} is generally similar to optimal target node set C * binary . This is likely due to the fact that the driver nodes which were randomly selected happened to coincide with the nodes which have high in-degrees, which results in reduced path distances between driver nodes and target nodes, and thus reduced control energy cost.
Finally, we apply TPGM to real networks and model networks of N=300, and compile the obtained results in Table 1, comparing with the control energy of initial random selection E (C 0 ) and random selection E (C rand ) . With the exception of Circuit-s838 and Teacher-student, where we drive the network with M/N = 0.5 and M/N = 0.6 fractional number of driver nodes to target control P/N = 0.7 and P/N = 0.75 fractional number of nodes, each of the network is driven by M/N = 0.4 fractional number of driver nodes to target control P/N = 0.6 fractional number of nodes. The energy cost of target node set selected by degree-based selection scheme is presented in Table 2.
The heuristic selection schemes lead to target control matrices C rand , C in.asc , C in.dsc , C out.asc , C out.dsc . Respectively, C rand corresponds to random selection of target nodes, where repeats are disallowed; C in.asc refers to target nodes chosen in ascending order according to their weighted in-degrees. Thus, when choosing target nodes for C in.asc , nodes with the lowest weighted in-degrees are selected. Likewise, C in.dsc is associated with target nodes chosen in descending order of weighted in-degrees, and nodes with the largest weighted in-degrees are selected. Analogously, C out.asc relates to target nodes picked in ascending order of weighted out-degrees; and C out.desc identifies target nodes of descending order of weighted out-degrees.
Examining Table 1, we observe that for any network, starting from its initial, E (C 0 ) , TPGM algorithm iteratively updates the C k matrix in the direction of quickest decreasing control energy until we obtain the convergent solution, E (C * ) , where the energy cost is typically reduced by a few orders of magnitude. Based on "Selecting binary optimal target node set from C * " section, we can choose the target nodes from the obtained optimal www.nature.com/scientificreports/ solution, C * , to obtain C * binary . While the conversion from dense real matrix solution to sparse binary matrix will in general result in increased control energy, comparing the energy cost E (C * binary ) to heuristic selection strategies of target nodes, such as E (C rand ) in Table 1 and degree-based selection in Table 2, we observe that the control energy compares favourably.

Discussion
When the learning rate η is chosen appropriately, TPGM is convergent. To be illustrative, the success rates of convergence as a function of η for different network topologies are shown in Fig. 5 Generally, when η is lower, the success rate is higher. Furthermore, when P is lower and the search space is smaller, higher η tend to be accommodated. All C * results presented in this work are convergent, even though the chosen η parameter may  Table 1. Control energy needed in various networks when different strategies are applied to selecting target control matrix C. E (C 0 ) is the average control energy from 10 independent initializations of random matrix without optimization. �E (C * )� is the average control energy from 10 independent iterative search using TPGM. We then convert these ten solutions into binary matrices E ([C * binary ] t ) , for t = 1, 2, . . . , 10 . The lowest of these is chosen as E (C * binary ) . E (C rand ) is the average control energy of 100 independent realizations of selecting target nodes randomly. Additional information such as standard deviations and mean of E ([C * binary ] t ) is presented in Supplementary Information. www.nature.com/scientificreports/ not have unity success rate. This is because, programmatically, for each independent iterative search, we can check the number of steps that the algorithm has taken; if the number of steps is abnormally small, the iteration is deemed to be non-convergent and to have terminated prematurely due to computation errors. Accordingly, we repeat the search iteration and in this way, convergence is always guaranteed. The convergence process is shown in Fig. 6a,b. In Fig. 6a, the k-th step cosine angle, cos(θ k ) , between energy gradient ∇E (C T k ) and projected energy gradient TCT k ∇E (C T k ) , approaches zero as TPGM iteration increases. The algorithm is deemed to have converged when the gradients are perpendicular. In Fig. 6b, we observe that each iteration of TPGM brings the matrix C k towards a lower energy cost, starting with a steep decrease in energy cost which becomes less steep with each iteration, until convergence.

Model
It is important to initialize C T 0 properly. In our experimentations, all elements of C T 0 are chosen randomly and set to be C T ij = 1 if node i is chosen to be the j-th target node initially, otherwise all other elements are set to be zero. We ensured row/column linear independence, which necessitates that no nodes are chosen repeatedly. If row/ column linear independence is not adhered, computations would not proceed due to mathematical error. While TPGM can accommodate initial C T 0 to be a random matrix with dense entries of random variables, this is Table 2. Control energy needed in various networks when degree-based selection strategies are applied to selecting target control matrix C. E (C in.asc ) refers to choosing target nodes in ascending order according to their weighted in-degrees. E (C in.dsc ) , E (C out.asc ) , and E (C out.dsc ) follows analogously, where dsc refers to order of descending and out refers to weighted out-degree.  Figure 5. The success rates of convergence as a function of log 10 η for ER10 network (of Fig. 2) for target controlling 3 and 5 nodes, and rhode food web network (of table 1) for target controlling P/N = 0.6 number of nodes. η varies from 1e−8 to 1e−3 . Each data point is calculated based on 100 independent iterative searches.
Scientific Reports | (2020) 10:18112 | https://doi.org/10.1038/s41598-020-75101-w www.nature.com/scientificreports/ ill-advised as the computation starting from a dense matrix is in general inefficient. Furthermore, the solutions C * and its associated C * binary also tend to perform poorly when compared to the former approach. However, no matter the initial condition, convergence is only affected by the learning rate parameter η.
When the number of nodes are small, such as in those shown in Figs. 1 and 2, we can systematically verify TPGM's searched solutions and compare them with the true optimal target node set, which can be found via a brute force exhaustive search approach. As N increases, the N choose P brute force search is no longer viable. However, such limitation should also be viewed as a strength, because it shows the difficulty of the search problem due to its infeasibly large search space; by relaxing matrix C into continuous search space, TPGM is able to search for C * in the direction of quickest energy cost decrease. Correspondingly, C * binary is recovered as the 'best signal' from C * .
The computational complexity of TPGM is O (N 3 ) owing to its iterative matrix multiplication. It should be noted that the computations of the controllability Gramian matrix, W = t f t 0 e Aτ BB T e A T τ dτ as well as the calling of the exponential function for computing e At f ( e A T t f ) are costly, which can result in bottlenecks when calculating them in the loop each time. An efficient way to code the TPGM is to compute the controllability Gramian W and e At f ( e A T t f ) outside the while loop and storing them as variables to be retrieved within the loop. Furthermore, instead of directly computing the integral, an efficient way to compute the controllability Gramian is to use the method of Ref. 41 .
The condition number of the controllability Gramian matrix plays an important role in our numerical experiments. When trying to target control a network using only driver nodes found from structural controllability 6 , the condition number of the Gramian may be too high and results in an infeasibly high control energy 31 . Increasing the number of control signals can lower the condition number and render the computation of the Gramian feasible 11,31 . Note that the condition number is dependent on matrices A, B, as well as time horizon t f . For some networks, despite increasing the number of driver nodes, the condition number of the Gramian may still be infeasibly high. To lower the condition number, we can normalise the link weights of the connection matrix, {a ij } , by dividing throughout a normalization constant. For networks whose connection strengths are not specified, we set {a ij } = 1 if there is a directed link from node j to node i. However, if this results in an infeasible condition number, then we will set the link strength to be random uniform [0.5, 1.5], which tends to make the condition number feasible.
Building a connection between target and driver node set optimization could be very interesting. Drawing on driver nodes optimization of previous works [15][16][17]22,30 for controlling the full node set, the conclusions formed are as follows: when driver nodes are placed in such a way such that they are evenly spaced out, the control energy is most optimal. Thus, by grouping controlled nodes to their nearest driver nodes and arranging the driver nodes in such a way that no groups have excessive geodesic paths 32 from drivers to controlled nodes, energy cost is most optimized. We expect the same conclusion to hold even for the case of target control, because full control and target control are not fundamentally different. Further, we expect that for the reverse problem of optimizing B, given that C is fixed, the same energy-optimal configurations, as presented in Figs. 1 and 2, to be found. In other words, whether optimizing driver nodes placement, or target nodes arrangement, the energy-optimal configurations are the same. Therefore, energy cost minimization could be an important mechanism in explaining structural self-organized configurations of driver/ target nodes in a natural or man-made complex system.
In conclusion, this paper demonstrate the possibility of minimizing control energy of directed complex networks by optimizing the target node matrix C, given that network connection matrix A, driver node placement B, and time horizon t f are fixed. We achieve this by adapting target control 19 , as well as deriving the energy gradient into the TPGM algorithm 22 . By ascribing target nodes to their nearest driver nodes, a directed complex network driven by driver nodes can be decomposed into elementary topologies 6,30 , Figure 6. An illustration of one particular iterative search's convergence process on the electronic circuit network, Circuit-s208. (a) cos(θ k ) moves closer to zero with iteration, although it is not always non-increasing. (b) Control energy is always non-increasing with each iteration, moving in the direction of largest decreasing energy cost. Insets show the first 0.3 × 10 5 iterations.

Scientific Reports
| (2020) 10:18112 | https://doi.org/10.1038/s41598-020-75101-w www.nature.com/scientificreports/ such as stem, circle, and dilation. Through extensive simulations on elementary topologies, our results reveal that control energy is most optimal when target nodes are chosen such that path distances from driver nodes to target nodes are minimized, corroborating existing literature results 21,33 . Furthermore, we validate our results on model networks (ER and SF) and real networks and show that optimal target node set (both C * and C * binary ) has control energy of a few orders of magnitude lower compared to target nodes chosen from heuristic selection schemes, such as random or nodes degree-based selection. Compared to previous works of control energy optimization 16,17,20,22,27 , which focus on driver node set optimization, our work is different in the sense that we instead considered target nodes selection optimization. Furthermore, our work show that in the context of target control, the choice of target nodes can account for variability of the control energy of a few orders of magnitude.
The problem of optimizing target nodes in the interest of control energy could be applicable to linear opinion networks [42][43][44] , where the state vector represents opinions of individuals, and the driver nodes are modeled as agents of influence. In such a system, we may be interested in influencing a certain fractional share of opinions to align with a pre-defined favorable opinion. Much like a voter model problem 45 , where we are only interested in obtaining the majority share of opinions, the specificity of which individuals to target control is not so much important as the control energy, which we would like to minimize. There are some avenues of research which appear to be promising. For example, within the framework of network controllability, recent works have incorporated conformity behavior 46,47 , where individuals' opinions adapt over time to mirror the average of their neighbors' , thus making the network dynamics richer and more realistic. It would be interesting to explore how optimal target nodes relate to a linear opinion network, both with and without conformity behavior in the future.

Methods
Model networks. Similar to recent works 18,19 , the model networks considered in this paper are modeled with stable dynamics, i.e. {a ii } < 0 ∀i . Specifically, we choose the diagonal of the connection matrix A to be chosen random uniformly from [−1, 1] , and then stabilize the nodal dynamics with {a ii } = δ i + ǫ , where ǫ is chosen such that the eigenvalues of A are all negative and the largest eigenvalue is −1 . The scale-free model network is constructed from the static model 48,49 , and the link weights {a ij } are drawn randomly from a uniform interval [0.5, 1.5].
Input nodes. The driver nodes selected are chosen using the method detailed in Ref. 6 , which uses the Hopcroft-Karp algorithm 7 to find the driver nodes to ensure controllability. To be specific, the driver nodes found in the research work presented in this paper are an overestimation, and guarantee full controllability of the complex network. Thus, controllability is always ensured, regardless of the choice of target node set C, and the term (CWC T ) is always invertible.
Weighted nodes degree. When considering nodes degree in the degree-based selection strategies, we computed the weighted link weights of each node to determine the in-degree (out-degree). For example, the indegree of node i is computed as: = N j=1,j� =i {a ij } . Note that self-links, {a ii } , do not count as node degree as they arise from categorically different sources 12 .
Gramian computation. The controllability Gramian can be efficiently calculated using the method of Ref. 41 : where each block partitioned matrix of (11) is a N × N matrix. The controllability Gramian is computed using (12).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.