Energy cost for controlling complex networks with linear dynamics

Examining the controllability of complex networks has received much attention recently. The focus of many studies is commonly trained on whether we can steer a system from an arbitrary initial state to any final state within finite time with admissible external inputs. In order to accomplish the control at the minimum cost, we must study how much control energy is needed to reach the desired state. At a given control distance between the initial and final states, existing results have offered the scaling behavior of lower bounds of the minimum energy in terms of the control time. However, to reach an arbitrary final state at a given control distance, the minimum energy is actually dominated by the upper bound, whose analytic expression still remains elusive. Here we theoretically show the scaling behavior of a precise upper bound of the minimum energy in terms of the time required to achieve control. Apart from validating the analytical results with numerical simulations, our findings are applicable to any number of nodes that receive inputs directly and any types of networks with linear dynamics. Moreover, more precise analytical results for the lower bound of the minimum energy are derived with the proposed method. Our results pave the way for implementing realistic control over various complex networks with the minimum control cost.

Figure 1

The lower and upper bounds of control energy for $n$ driver nodes. By controlling all nodes directly, here we show the numerical and analytical results for lower ($\underline{E}$) and upper ($\overline{E}$) bounds of control energy for different types of $\mathbf{A}$. To adjust the maximum (minimum) eigenvalue of $\mathbf{A}$ intuitively, we set the link weight $a_{ij}$ uniformly from $[0,1]$ in (a) to (d) and from $[-1,0]$ in (e) and (f); each self-loop (diagonal element) is set as $a+s_i$ with $s_i=-{\sum }_{j=1}^n{a}_{ij}$. In (a), we set $a=-5$, which guarantees $\mathbf{A}$ is ND with eigenvalues in $[-14.0266,-5]$. Similarly, in (b), $a=0$ and $\mathbf{A}$ is NSD with eigenvalues in $[-8.5243,0]$. In (c) and (d), we have $a=5$, and $\mathbf{A}$ is ID with eigenvalues in $[-4.0266,5]$. In (e), we set $a=0$, and hence $\mathbf{A}$ is PSD with all eigenvalues in $[0,8.3062]$. In (f), $a=5$ and $\mathbf{A}$ is PD with all eigenvalues in $[5,13.7144]$. In each panel, triangles (blue and purple) represent results obtained by numerical calculations and full lines indicate analytical derivations under our framework (see Sec. 3a and Table 1). For small $t_f$, from each panel with horizontal axis $\text{ln}\left(t_f\right)$, we see that all slopes are $-1$, which confirm our analytical results that both $\overline{E}$ and $\underline{E}$ approximate $\frac{1}{t_f}$ for different types of $\mathbf{A}$. For large $t_f$, subgraphs with horizontal axis $t_f$ or ln($t_f$) show the analytical scaling behaviors of the bounds of energy precisely. Here we adopt the BA scale-free network with $n=50$, and network is constructed based on the preferential attachment with average degree 5.8 [27].

Figure 2

Veracity of eigenvalues estimation based on Eqs. (5) and (6). In (a), we randomly generate 25 matrices with minimum eigenvalue being $4i, i=1,2,\cdots ,25$, where $i$ is the index of the matrix. The horizontal and vertical coordinates represent the true eigenvalues and estimated eigenvalues by Eqs. (5) and (6), from which it is clear the generated pattern almost overlaps with $y=x$. The inset presents ratio errors of differences between approximated eigenvalues by Eqs. (5) and (6) and the true eigenvalues, which indicates the accuracy of estimation is reliable, especially the estimation of minimum eigenvalues by (6). In (b) and (c), we make estimations of the maximum and the minimum eigenvalues of the matrix $\mathbf{G}$ for the case of $d$ driver nodes. All validations are carried out on BA networks with 10 nodes. To be more persuasive, we set the number of driver nodes denoted by $d$ as $1,2,4,6,8$, respectively. And for BA networks, all cases of the different $\mathbf{A}$ with different properties (ND, NSD, ID, PD, PSD) are considered. In (b), the horizontal axis represents the real minimum eigenvalue of $\mathbf{M}$, and the vertical axis represents the estimated value by (6), where we set $\mathbf{A}$ as PD, PSD, not PD by selecting $a_{ij}$ from $[-4,-1]$ uniformly and $a$ from $-4$ to 4 with an interval 0.2. In (c), the horizontal axis represents the real maximum eigenvalue of $\mathbf{M}$, and the vertical axis represents the estimated value by (5), where we set $\mathbf{A}$ as ND, NSD, not ND by uniformly selecting $a_{ij}$ from $[1,3]$ and $a$ from $-4$ to 4 with an interval 0.2. Subplots present the ratio error similar to those in (a) with half original data, where the horizontal axis presents the matrix index and the vertical axis indicates the ratio error.

Figure 3

The lower bound of energy comparisons between the methods shown in Ref. [16] and this paper. Here we randomly generate BA scale-free networks with $\mathbf{A}$ being ND (other parameters are the same as those in Fig. 1) and $a_{ij}$ is selected from $[1,3]$ uniformly with $a=-2$. For approximating the maximum eigenvalue of $\mathbf{M}$, here we use the method shown in (5), while in Ref. [16], it is inferred by the corresponding trace. Since the existing results only consider the scenario for one driver node, we follow this setting. In (a), the network size is set as 10. In (b), the analytically derived and numerical lower bounds of $E$ are depicted at $t_f=200$ for different network sizes ($n$) accordingly. For all cases, we can see that the method we employed generates much more precise $\underline{E}$ compared to the existed tools.

Figure 4

The lower and upper bounds of energy for one driver node. The scaling behavior of the lower and upper bounds of energy cost is given for one driver node, and the summation of analytical results are presented in Tables 1 and 2. For scale-free networks in (a)–(f) and random networks in (g)–(l), all parameters are the same except that the network structure is different. In (a)–(c), with small $t_f, \underline{E}\sim t_f^{-1}$ for all $\mathbf{A}$. In (d)–(f) for upper bound, the slope of triangular trajectory is much less than $-1$. Parameters are selected the same as those given in Fig. 1. The interval of the uniform distribution is $[0,1]$ in (a)–(c), $[1,3]$ in (d), $[-1,0]$ in (e), and $[-5,-2]$ in (f). In (a), $a=-5$, by which $\mathbf{A}$ is ND with eigenvalues in $[-14.0266,-5]$. Similarly, in (b) and (e), $a=0$ such that $\mathbf{A}$ is NSD and PSD, respectively. In (c) and (d), $a=5$ such that $\mathbf{A}$ is ID. In (f), $a=3$, such that the minimum eigenvalue of $\mathbf{A}$ is 3. In (g)–(l), the probability for adding an edge between every randomly selected pair nodes is 0.1 [29].

Figure 5

The lower and upper bounds of control energy for 20 driver nodes. In (a)–(c), with small $t_f, \underline{E}\sim t_f^{-1}$ for all $\mathbf{A}$. In (d)–(f) for upper bound, the slope of triangular trajectory is much less than $-1$. The summation of the analytical results are presented in Tables 1 and 2. Parameters are selected as those given in Fig. 1. The interval of uniform distribution is $[0,1]$ in (a)–(d), and $[-1,0]$ in (e) and (f). In (a), $a=-5$, by which $\mathbf{A}$ is ND with eigenvalues in $[-12.5048,-5]$. Similarly, in (b) and (e), $a=0$ such that $\mathbf{A}$ is NSD and PSD, respectively. In (c) and (d), $a=5$ such that $\mathbf{A}$ is ID. Similarly, $a=5$ such that $\mathbf{A}$ is PD.

Figure 6

The upper bound of energy estimated by (7) with ${\mathbf{M}}_1$ representing $\mathbf{M}$. All networks are BA networks with 50 nodes and one driver node. In (a) and (b), with large $t_f, \overline{E}$ converges to the constant when $\mathbf{A}$ is NSD or ID. And the trajectories of $\overline{E}$ with $\mathbf{M}$ and ${\mathbf{M}}_1$ handled are almost coincide. Parameters are selected as those given in Fig. 1. In (a), we set $\mathbf{A}$ as NSD by setting $a=0$, where the eigenvalues of $\mathbf{A}$ locate in $[-9.0266,0]$. In (b), we generate BA network with $\mathbf{A}$ being ID by setting $a=5$, where the eigenvalues of $\mathbf{A}$ locate in $[-5.0028,5]$.