Minimum energy control for complex networks

The aim of this paper is to shed light on the problem of controlling a complex network with minimal control energy. We show first that the control energy depends on the time constant of the modes of the network, and that the closer the eigenvalues are to the imaginary axis of the complex plane, the less energy is required for complete controllability. In the limit case of networks having all purely imaginary eigenvalues (e.g. networks of coupled harmonic oscillators), several constructive algorithms for minimum control energy driver node selection are developed. A general heuristic principle valid for any directed network is also proposed: the overall cost of controlling a network is reduced when the controls are concentrated on the nodes with highest ratio of weighted outdegree vs indegree.

where x ∈ R n is the state vector, A ∈ R n×n is the state update matrix, B ∈ R n×m is the input matrix, and u is the m-dimensional input vector. The reachable set of (S1) in time t f from x o is the set where φ(t, u, x o ) is the solution of (S1) at time t with input u and Ω is the admissible set of control inputs, here Ω = R m . The system (S1) is reachable (or controllable from the origin [1]) in time t f if any x f ∈ R n can be reached from 0 by some control u ∈ Ω in time t f , i.e. if R t f (0) = R n . It is said controllable to the origin if any x o ∈ R n can be brought to 0 by some control u ∈ Ω in time t f . The system (S1) is said Finite-time Gramians. The time-t f reachability (or controllability from 0) Gramian is the symmetric matrix while the time-t f controllability to 0 Gramian (normally called the controllability Gramian, [1]) is The two Gramians are positive definite whenever (A, B) is controllable, and are related by Similarly, for their inverses, Finite-time control energy for state transfer. The transfer of the state from any x o into any other x f in time t f can be accomplished by many controls. In order to quantify how costly a state transfer is on a system, one can choose to consider the control input that minimizes the input energy, i.e., the functional Such control can be computed explicitly [1] as and the corresponding transfer cost as To be more precise on how the control energy is formed, we have to split the state transfer energy (S6) into subtasks: 2. x f = 0 (controllability to 0 problem) In particular, both W r (t f ) and W c (t f ) enter into the cost function. Notice that, to quantify the amount of control energy of these problems we need to compute the inverse of W r (t f ) and W c (t f ). Let us look at how the stability/instability of the eigenvalues influences the two costs E r (t f ) and E c (t f ).
• If A is stable (i.e., Re[λ(A)] < 0), then "escaping from 0" (i.e., the reachability problem) requires more energy than transferring to 0, (i.e., the controllability to 0 problem) because the modes of A naturally tend to converge to 0.
• If A is antistable (i.e., Re[λ(A)] > 0, −A is stable), then the opposite considerations are valid: the modes of A tend to amplify the magnitude of the state, simplifying the reachability problem but complicating the controllability to 0 problem.
• If A has eigenvalues with both negative and positive real part, the two situations coexist.
Hence computing a plausible measure of control energy for a generic state transfer x o → x f requires to take into account the "difficult" directions of both cases.
Gramian-based measures of control energy. Expression like E(t f ) depend on the choice of x o and x f . In control theory, it has been known for a long time that it is possible to obtain more objective estimates of the control energy by computing figures of merits based on the Gramian [18]. The following metrics are often used: 1. λ min (W r ) = 1 λmax(W −1 r ) : the min eigenvalue of the Gramian is a worstcase metric, equal to the inverse of the max eigenvalue of W −1 r . λ max (W −1 r ) estimates the energy required to move along the direction which is most difficult to control, hence when λ max (W −1 r ) descreases the control energy decreases. In turn, the control energy decreases when λ min (W r ) increases.
2. tr(W r ): the trace of the Gramian is inversely proportional to the average energy required to control a system, hence when tr(W r ) increases the control energy decreases.
3. tr(W −1 r ): the trace of the inverse of the Gramian is proportional to the average energy needed to control the system, i.e., the control energy decreases with tr(W −1 r ).
In summary, minimizing the control energy means maximizing the first and second measure or minimizing the third.

Control energy: infinite-time horizon formulation
When t f → ∞, then E(t f ) converges (or diverges) to a quantity and so do E r (t f ) and E c (t f ).
When t f → ∞, both Gramians become infinite-time integrals, which may be convergent or divergent, depending on the modes of A. If A stable, then to exist finite and positive definite when (A, B) controllable. In the mixed eigenvalues cases the two expression (S8) and (S9) both diverge, although their inverses both exist (see below): W −1 r depends only on the stable modes, and W −1 c only on the unstable modes. In an expression like (S7), for t f → ∞ the "cross-terms" . They also vanish when A is antistable because the same cross terms can also be written as Controllability to 0 in the infinite-time horizon. Let us observe what happens for instance to the controllability to 0 problem according to the eigenvalues of A.
• If A is stable, then in correspondence of u = 0, lim t→∞ x(t) = 0 for all x o , meaning that the controllability to 0 problem can be solved with zero energy E c = 0. Furthermore, the integral (S2) converges to (S8), whose value can be computed solving the following Lyapunov equation: Such a solution always exists and it is W r > 0 (positive definite) if the pair (A, B) is controllable.
• When instead A is antistable, then the integral (S2) diverges as t f → ∞. Hence the solution with u = 0 is no longer feasible as all modes are unstable (and diverge as soon as x o = 0), meaning that to find a minimizer of (S7) we have to proceed in some other way. Since (A, B) controllable, we can determine u(t) as if we were computing a stabilizing feedback law, i.e., expressing u(t) as a function of the state x(t) so that the resulting closed loop system converges to 0 asymptotically. Such a feedback law can be computed solving in P an algebraic Riccati equation (ARE) Such ARE admits a positive definite solution P , which can in turn be computed solving in L the Lyapunov equation (in −A, which is stable, hence a solution L > 0 always exists) and then setting P = L −1 . It can be verified directly that the controllability Gramian W c in (S9) is one such solution L, i.e., L = W c solves (S12). Correspondingly we obtain P = W −1 c . From the theory of linear-quadratic regulators (in particular [24], Ch. 10) the feedback controller guarantees stability of the closed-loop systeṁ i.e., A − BB T P is a stable matrix. The feedback law (S13) also minimizes the input energy (S7) which is equal to • When A has eigenvalues with both positive and negative real part and no purely imaginary eigenvalues, then the two situations described above occur simultaneously. Assume A is split into two diagonal blocks, one consisting of only eigenvalues with negative real part and the second only of eigenvalues of positive real part. This can always be achieved through a change of basis [27]. Split B and x(t) accordingly: > 0 (S14) In the infinite time horizon, the u = 0 control steers optimally the x 1 subvector, while for the x 2 part a feedback controller provides the energy-minimizing solution. From (S11) we obtain that the ARE has solution where P 2 solves the ARE for the (A 2 , B 2 ) subsystem. Hence the control input achieves a transfer to the origin with minimal energy cost equal to Furthermore, combining (S10) and (S12), we have that for (S14) the following two Lyapunov equations must hold simultaneously: Reachability in the infinite-time horizon. Let us now consider the reachability problem (i.e., controllability from 0). Now the roles of stable and unstable eigenvalues are exchanged.
• If A stable, reachability requires an active control (here a destabilizing state feedback law) in order to steer x(t) out of the origin. The energy-optimal solution consists in choosing u = −B T P x(t) with P > 0 solution of the ARE P A + A T P + P BB T P = 0 (S16) or, equivalently, P = K −1 with K solution of the Lyapunov equation A is stable, hence K > 0 solving (S17) and P > 0 solving (S16) always exist. The resulting closed loop matrix A − BB T P must be antistable. From (S10) and (S17) it can also be K = W r , the reachability Gramian.
• If A antistable, then u = 0 is the minimal energy controller (an infinitesimal amount energy at t = 0 is enough to "kick" the system towards the right direction x f when initialized in x o = 0; this amount of energy is negligible in the infinite time horizon considered here). Since −A stable, a Lyapunov equation like (S12) holds, with solution L = W c .
• When A has eigenvalues with both positive and negative real part and no purely imaginary eigenvalues, then a decomposition like (S14) can be obtained through a change of basis. The complete ARE has now solution and the controller achieving the transfer with minimal energy is for an amount of energy equal to The decomposition (S14) also in this case induce a pair of Lyapunov equations identical to (S15).

Mixed Gramian in infinite and finite time horizon
In order to assemble the considerations of the previous sections, it is useful to introduce a third Gramian which we call mixed Gramian, W m and which gathers the directions difficult to control of both the reachability and the controllability to 0 problems. This mixed Gramian is used extensively in the paper.
Infinite-time horizon mixed Gramian Assume that the spectrum of A contains k, 0 k n, eigenvalues with negative real part, and n − k eigenvalues with positive real part (and no purely imaginary eigenvalues). Then, as already mentioned above, there exist a change of basis V bringing A into the form (S14): and, correspondingly, with Re[λ(Ā 1 )] < 0 and Re[λ(Ā 2 )] > 0. In the new basis, the two Lyapunov equations (S15) hold, which can be rewritten as the mixed Gramian. Following [27], the expression of the mixed Gramian in the original basis is By construction, the mixed Gramian matrix W m always exists when A has no purely imaginary eigenvalues, and it summarizes the infinite-horizon contribution of the stable eigenvalues to the reachability problem and of the unstable eigenvalues to the controllability to 0 problem (i.e., all cases leading to a high control energy).
Control energy measures for the mixed Gramian The three control energy metrics introduced above can be computed also for the mixed Gramian, with an analogous interpretation. In fact λ min (W m ), tr(W m ) and tr(W −1 m ) are the three measures that are mostly used in the paper to quantify the control energy.
Finite-time horizon mixed Gramian Using the insight given by the previous arguments, it is possible to construct also a finite-time mixed Gramian, which weights only the modes that are difficult to control in the two state transfer problems. In the basis in which A is split into stable and antistable diagonal blocks, (S18), this is given bȳ are the equivalent of (S2) and (S3) for the two subsystems (Ā 1 ,B 1 ) and (Ā 2 ,B 2 ). An equation like (S19) has no coupling terms between the two subsystems (i.e., terms of the formB 1B2 ). These terms disappear asymptotically, but transiently they give a contribution, hence the finite-time formulation ofW m (t f ) is only an approximation. In the original basis,

and the input energy is
Clearly a proxy for this quantity is obtained by simply flipping the sign the real part of the unstable eigenvalues of A and considering only the reachability problem on the resulting stable system. Equivalently, all eigenvalues can be made unstable, and the controllability to 0 problem considered. More details and explicit formulas for the finite-time horizon are provided in e.g. [14,22].

Controllability with bounded controls
Consider the system (S1). Assume u ∈ Ω, where Ω is a compact set of R m containing the origin in its interior. Assume (A, B) is controllable. Then we have the following, see [13,5] and [23], p. 122.
• A necessary and sufficient condition for the origin to be steered to any point of R n in finite time (i.e., the reachability problem) is that no eigenvalue of A has negative real part.
• A necessary and sufficient condition for any point of R n to be steered to the origin in finite time (i.e., the controllability to 0 problem) is that no eigenvalue of A has positive real part.
Combining the two: • A necessary and sufficient condition for complete controllability (from any point x o to any point x f ) in finite time is that all eigenvalues have zero real part.

Control of coupled harmonic oscillators
A network of n coupled harmonic oscillators can be written as a system of second order differential equations where M ii > 0 is the mass of the i-th oscillator, k i 0 its stiffness, k ij 0 the coupling stiffness between the i-th and j-th oscillators, and β i ∈ {0, 1} indicates the presence or absence of a forcing term in the i-th oscillator. In matrix form, (S20) can be rewritten as where M = M T = diag(M ii ) > 0 is the mass matrix, K = K T 0 the stiffness matrix, and B is a n × m matrix whose columns are the elementary vectors corresponding to the β i = 1. When u = 0, the solutions of (S21) have the form q = φe iωt in correspondence of the pairs ω j and φ j , j = 1, . . . , n, that are the solutions of the generalized eigenvalues/eigenvector equation The ω j are called the natural frequencies of (S21). Denote the matrix of eigenvectors. Φ can be used to pass to a so-called modal basis, in which the oscillators are decoupled. In fact, it can be verified directly that in correspondence of the change of basis has decoupled dynamics (but coupled inputs). The state space representation of (S21) is 2n dimensional. If In terms of the state space model (S24), the eigenvalues are λ j = ±iω j , j = 1, . . . , n, of eigenvectors If Ψ = M Φ, the state space representation in the modal basis A key advantage of the modal representation is that the reachability Gramian of the pair (A 1 , B 1 ) can be computed explicitly. As a matter of fact, when the eigenvalues are on the imaginary axis, there is no need to resort to the mixed Gramian introduced earlier. As the integral (S2) (or (S3)) diverges, the infinite-time Gramian cannot be computed. However, in the modal basis z, W z (t f ) (or more precisely W z,r (t f )) is diagonally dominant, and for t f sufficiently long it can be approximated by its diagonal terms. These terms are computed explicitly in [2]: If we assume that the mass matrix M is diagonal, then it is always possible to choose Ψ so that M 1 = I and K 1 diagonal, by suitably rescaling the eigenvectors ψ j . In this case and the Gramian is determined by the lower diagonal block. When the columns of B are elementary vectors as in our case, the product Ψ T M −1 BB T M −1 Ψ can be written explicitly as sum of rank-1 matrices: (only m of the n factors β j ∈ {0, 1} are nonzero) and its diagonal entries are Hence the expression for the Gramian in the modal basis is Notice the linearity in t f , meaning that all components diverge to ∞ with the same speed when t f → ∞. This expression can be used to compute the various measures of control energy we have adopted in the paper, and hence to optimize the driver node placement problem. For instance, selecting inputs according to λ min (W z ) amounts to solving the following MILP maxmin problem: which can be solved exactly only for systems of moderate size. However, efficient heuristics can be derived for it, such as Algorithm 1.
Algorithm 1 Driver node placement that maximizes λ min (W z ). Input: • y s = yˆi If instead we choose to maximize the tr(W z ), then we get Since tr(W z ) is linear in the β j , this is a linear optimization problem, hence solvable exactly and efficiently for any n. Finally, also for the minimization an efficient heuristic can be set, as outlined in Algorithm 2.
Algorithm 2 Driver node placement that minimizes tr(W −1 z ). Input: Looking at an expression like (S26), it is possible to understand what kind of behavior yields good controllability properties to certain driver nodes. For instance when measuring according to λ min (W z ), from (S26), the columns of Ψ T express how nodes in the original basis are spread among the state variables in the modal basis. What is needed is an eigenbasis such that the j-th component of all eigenvectors has "support" on all the directions of the state space, i.e., all ψ 1 j , . . . , ψ n j are nonvanishing and possibly all as large as possible. Since W z is approximated well by a diagonal matrix, the "coverage" effect of control inputs is additive, hence choosing a pair of controls i and j for which the sum of {ψ k i } k=1...,n and {ψ k j } k=1...,n has all components that are as large as possible guarantees an improvement in the control cost with respect to taking only one of the two inputs.
Notice that although M and K are both symmetric, KM −1 need not be, hence Ψ need not be an orthogonal matrix. It can be rendered orthogonal if a slightly different modal basis is chosen, see [10] for the details. In that case Ψ T = Ψ −1 i.e., the "coverage" discussed here is given by the left eigenvectors of (S22), condition sometimes considered in the literature [11].
Another basis for the state space that can be used in place of (S23) is given byx = qq . With this choice, the state space realization iṡ It is straightforward to verify that in this basis the roles of w in and w out are exchanged, hence a criterion for driver node selection becomes ranking according to w in /w out (instead of w out /w in ).

Datasets
The majority of the networks listed in Table 1 of the paper is composed of a large connected component, plus a number of other very small connected components. In our study, only the largest connected component is considered. The following networks are considered: • Biology, transcriptional networks -E.coli-transcr.: Gene regulatory network of E.coli, downloaded from RegulonDB database (http://regulondb.ccg.unam.mx). See [9].
-Yeast-metab: Metabolic network of S.cerevisiae. Assembled from the list of reactions in [8].
• Biology, signaling networks -Macrophage: The molecular interaction map of a macrophage obtained by [19].
-Toll-like: Signaling network for the Toll-like receptor. Assembled from [20].
• Eoclogical, Food-webs -Foodweb-Florida: Trophic dynamics in a south Florida ecosystem, from Pajek collection, see also [25]. -Advogato: Advogato is an online community platform for developers of free software launched in 1999. Nodes are users of Advogato and the directed edges represent trust relationships, see [15].
• Transport and Trade -US Airport: This is the directed network of flights between US airports in 2010. Each edge represents a connection from one airport to another. Compiled from the blog post "Why Anchorage is not (that) important: Binary ties and Sample selection", and downloadable from https://toreopsahl.com/datasets/#usairports. The power networks used in the second part of the paper are listed in Table S1. All nodes are treated equally, regardless of their function as generators or loads in the real grid.      Figure S3: Driver node placement for ER networks with p = 0.05. (a): Comparison between the value of the various measures of control energy obtained for driver node placement strategies based on r w = w out /w in (red, labelled "ordered") and the same measure for random driver node assignments (blue, labelled "random"). As can be seen on the ratios shown in (b), all measures improve, especially when m is low. Notice that also the condition number of W m (i.e., λ max (W m )/λ min (W m ), lower right panel) improves.  Figure S5: Driver node placement for SF networks with γ in = 3.14 and γ out = 2.87.
(a): Comparison between the value of the various measures of control energy obtained for driver node placement strategies based on r w = w out /w in (red, labelled "ordered") and the same measure for random driver node assignments (blue, labelled "random"). As the ratios in (b) show, the improvement in all measures is normally of several orders of magnitude. Also the condition number of W m (i.e., λ max (W m )/λ min (W m )) improves substantially.    . Red: driver node placement based on λ min (W z ). Violet: placement based on tr(W z ). Green: placement based on tr(W −1 z ). Cyan: placement based on w out /w in . Blue: random input assignment. All driver node placement strategies still beat a random assignment, but with worse performances with respect to Fig. 4. Of the four measures, λ min (W z ) and tr(W −1 z ) tend to behave similarly and so do w out /w in and tr(W z ) (in the mid plot they completely overlap, and both give the true optimum). (b): Overlap in the node ranking of the different driver node placement strategies. Color code is the same as in (a). The only highly significant overlap is still between w out /w in and tr(W z ) (> 90%), while λ min (W z ) and tr(W −1 z ) correspond to different node ranking patterns. None of the strategies orders nodes according to w in /w out , as expected.