arXiv:1612.03082v1 [] 9 Dec 2016

Minimum energy control for complex
networks
Gustav Lindmark and Claudio Altafini∗
Division of Automatic Control, Dept. of Electrical Engineering,
Linköping University, SE-58183, Linköping, Sweden.
March 12, 2018

Abstract
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.

Significance statement Controlling a complex network, i.e., steering the
state variables associated to the nodes of the networks from a configuration
to another, can cost a lot of energy. The problem studied in the paper
is how to choose the controls so as to minimize the overall cost of a state
transfer. It turns out that the optimal strategy for minimum energy control
of a complex network consists in placing the control inputs on the nodes that
have the highest skewness in their degree distributions, i.e., the highest ratio
between their weighted outdegree and indegree.
∗

Corresponding author: C. Altafini. Email: claudio.altafini@liu.se

1

1

Introduction

Understanding the basic principles that allow to control a complex network is
a key prerequisite in order to move from a passive observation of its functioning to the active enforcement of a desired behavior. Such an understanding
has grown considerably in recent years. For instance the authors of [22] have
used classical control-theoretical notions like structural controllability to determine a minimal number of driver nodes (i.e., nodes of the network which
must be endowed with control authority) that guarantee controllability of a
network. Several works have explored the topological properties underlying
such notions of controllability [8, 11, 21, 26, 27, 31], or have suggested to
use other alternative controllability conditions [10, 25, 47]. Several of these
approaches are constructive, in the sense that they provide receipts on how
to identify a subset of driver nodes that guarantees controllability. However, as observed for instance in [45], controllability is intrinsically a yes/no
concept that does not take into account the effort needed to control a network. A consequence is that even if a network is controllable with a certain
set of driver nodes, the control energy that those nodes require may result
unrealistically large. Achieving “controllability in practice” i.e., with a limited control effort, is a more difficult task, little understood in terms of the
underlying system dynamics of a network. In addition, in spite of the numerous attempts [4, 7, 20, 25, 28, 29, 37, 38, 41, 45, 46], no clear strategy
has yet emerged for the related problem of selecting the driver nodes so as
to minimize the control energy.
The aim of this paper is to tackle exactly these two issues, namely: i)
to shed light on what are the dynamical properties of a network that make
its controllability costly; and ii) to develop driver node placement strategies
requiring minimum control energy. We show in the paper that for linear dynamics the natural time constants of the modes of the system are key factors
in determining how much energy a control must use. Since the time constants
of a linear system are inversely proportional to the real part of its eigenvalues,
systems that have eigenvalues near the imaginary axis (i.e., nearly oscillatory behavior) are easier to control than systems having fast dynamics (i.e.,
eigenvalues with large real part). For networks of coupled harmonic oscillators, which have purely imaginary eigenvalues, we show that it is possible to
obtain explicit criteria for minimum energy driver node placement. One of
these criteria ranks nodes according to the ratio between weighted outdegree
and weighted indegree. We show that for any given network such criterion
2

systematically outperforms a random driver node assignment even by orders
of magnitude, regardless of the metric used to quantify the control energy.

2

Methods

Reachability v.s. Controllability to 0. A linear system
ẋ = Ax + Bu

(1)

is controllable if there exists an input u(t) that transfers the n-dimensional
state vector x(t) from any point xo to any other point xf in Rn . The Kalman
rank condition for controllability, rank([B AB A2 B . . . Ak B]) = n for k sufficiently large, only provides a yes/no answer but does not quantifies what
is the cost, in term of input effort, of such state transfer. A possible approach to investigate “controllability in practice” consists in quantifying the
least energy that a control requires to accomplish the state transfer, i.e.,
in computing u(t) mapping xo in xf in a certain time tf while minimizing
Rt
E(tf ) = 0 f ku(τ )k2 dτ . For linear systems like (1), a closed form solution to
this problem exists and the resulting cost is
E(tf ) = (xf − eAtf xo )T Wr−1 (tf )(xf − eAtf xo ),

(2)

Rt
T
where the matrix Wr (tf ) = 0 f eAτ BB T eA τ dτ is called the reachability
Gramian [2]. The control that achieves the state transfer xo → xf with
minimal cost can be computed explicitly:
u(t) = B T eA

T (t −t)
f

Wr−1 (tf )(xf − E Atf xo ),

t ∈ [0, tf ].

(3)

Various metrics have been proposed to quantify the difficulty of the state
transfer based on the Gramian, like its minimum eigenvalue λmin (Wr ), its
trace tr(Wr ), or the trace of its inverse tr(Wr−1 ), see [24] and also SI for a
more detailed description.
We would like now to describe how (2) depends on the eigenvalues of
A. In order to do that, one must observe that (2) is the sum of contributions originating from two distinct problems: 1): controllablity-from0 (or reachability, as it is normally called in control theory [2]) and 2):
controllablity-to-0. The first problem consists in choosing xo = 0, in which
case (2) reduces to Er (tf ) = xTf Wr−1 (tf )xf , while in the second xf = 0 leads
3

T

to Ec (tf ) = xTo Wc−1 (tf )xo where Wc (tf ) = e−A tf Wr (tf )e−Atf is a second
Gramian, called the controllability Gramian. The two problems are characterized by different types of difficulties when doing a state transfer, all
related to the stability of the eigenvalues of A. For instance the reachability
problem is difficult along the stable eigendirections of A because the control
has to win the natural decay of the unforced system to 0, while the unstable
eigenvalues help the system escaping from 0 by amplifying any small input
on the unstable eigenspaces, see Fig. 1 for a graphical explanation. The surfaces of Er (tf ) shown in Fig. 1 (a) reflect these qualitative differences. On
the contrary, the influence of the eigenvalues of A is the opposite for the
controllability-to-0 problem shown in Fig. 1 (b). Hence if we want to evaluate the worst-case cost of a transfer between any xo and any xf (problem
sometimes referred to as complete controllability [35]), we have to combine
the difficult cases of the two situations just described. This can be done
combining the two Gramians into a “mixed” Gramian Wm obtained splitting
A into its stable and antistable parts and forming a reachability subGramian
for the former and a controllability subGramian for the latter, see SI for the
details. Such Gramian can be computed in closed form only when the time
of the transfer tends to infinity. In the infinite time horizon, in fact, both
Wr and Wc diverge, but their inverses are well-posed and depend only on the
stable modes the former and the unstable modes the latter. These are the
parts constituting the inverse of Wm , see Fig. 1 (c). A finite-horizon version
of Wm (and Wm−1 ) can be deduced from the infinite horizon ones.
Eigenvalues of random matrices.
The so-called circular law states that
√
a matrix A of entries aij / n where aij are i.i.d. random variables with
zero-mean and unit variance has spectral distribution which converges to the
uniform distribution on the unit disk as n → ∞, regardless of the probability
distribution from which the aij are drawn [1]. A numerical example is shown
in Fig. 2(a) (top left panel). By suitably changing the diagonal entries, the
unit disk can be shifted horizontally at will, for instance rendering the entire
spectrum stable (Fig. 2(a), top mid panel) or antistable (Fig. 2(a), top right
panel). A random matrix is typically a full matrix, meaning that the underlying graph of interactions is fully connected. The circular law is however valid
also for sparse matrices, for instance for Erdős-Rényi (ER) topologies. If p is
√
the edge probability, then A = (aij )/ p · n still has eigenvalues distributed
uniformly on the unit disk, see Fig. S1(a).

4

A generalization of the circular law is the elliptic law, in which the unit
disk containing the eigenvalues of A is squeezed in one of the two axes. To
do so, the pairs of entries {aij , aji } of A have to be drawn from a bivariate
distribution with zero marginal means and covariance matrix expressing the
compression of one of the two axes, see [1]. Various examples of elliptic laws
are shown in the lower panels of Fig. 2 (a). Also elliptic laws generalize to
sparse matrices, see Fig. S1(a).

3

Results

Control energy as a function of the real part of the eigenvalues of A.
In a driver node placement problem, the inputs affect a single node, hence
the columns of B are elementary vectors, i.e., vectors having one entry equal
to 1 and the rest equal to 0. When A is a random matrix, the underlying
graph is generically fully connected, hence issues like selection of the number
of driver nodes based on the connectivity of the graph become irrelevant.
Having disentangled the problem from topological aspects, the dependence
of the control effort from other factors, like the spectrum of A, becomes
more evident and easier to investigate. If for instance we place driver nodes at
random and use the mixed Gramian Wm to form the various energy measures
mentioned above for quantifying the control effort, then we have the results
shown in Fig. 2(b). As expected, all indicators improve with the number of
inputs. What is more interesting is that when we repeat the computation
for the various spectral distributions of Fig. 2(a), the result is that the cost
of controllability decreases when the (absolute value of the) real part of the
eigenvalues of A decreases. All measures give an unanimous answer on this
dependence, regardless of the number of inputs considered. In particular,
when A has eigenvalues which are very close to the imaginary axis (lower
right panel of Fig. 2(a) and cyan curves in Fig. 2(b)) then the worst-case
controllability direction is easiest to control (i.e., λmin (Wm ) is bigger), but
also the average energy needed for controllability on all directions decreases
(i.e., tr(Wm ) increases and tr(Wm−1 ) decreases).
Recall that in a linear unforced dynamical system the real part of the
eigenvalues of A determines how fast/slow a system converges to the origin (stable eigenvalues, when real part of λ(A) is negative) or diverges to
∞ (unstable eigenvalues, when real part of λ(A) is positive). Such convergence/divergence speed grows with the absolute value of the real part of λ(A).
5

In the complete controllability problem, both stable and unstable modes of
A are gathered together, and all have to be “dominated” by the controls to
achieve controllability. When the modes of the system are all slow, like when
they are very close to the imaginary axis, then the energy needed to dominate them all is lower than when some of them are fast (i.e., the eigenvalues
have large real part).
An identical result is valid also for sparse matrices. In particular, for ER
graphs with edge probability p = 0.05 and coefficients from a bivariate normal
distribution (yielding elliptic laws as in Fig. S1(a)), the various norms used
to quantify input energy are shown in Fig. S1(b). Their pattern is identical
to the full graph case of Fig. 2(b).
The computations shown in Fig. 2(b) are performed with the infinitehorizon mixed Gramian Wm described in the SI, because such Wm can be
easily computed in closed form. A finite-horizon Wm (tf ) can be approximately obtained from it, but the arbitrarity of tf makes it hard to set up an
unbiased comparison of the various spectral distributions of A of Fig. 2(a),
which are characterized by widely different time constants (inversely correlated to the amplitude of the real part of λ(A)). Observe in Fig. S2 how the
various measures of controllability computed with a finite-time Wm (tf ) tend
all to the infinite-time Wm but with different speeds.
Driver node placement based on weighted connectivity. In the analysis carried out so far the driver nodes were chosen randomly. A topic that
has raised a remarkable interest in recent times (and which is still open in
our knowledge) is devising driver node placement strategies that are efficient
in terms of input energy [4, 7, 20, 28, 29, 37, 41, 46]. If we consider as
weighted indegree and outdegree at node i the sum of the weights
in abPn
solute value
Pnof all incoming or outgoing edges, i.e.,√win (i) = j=1 |aij | and
wout (i) = j=1 |aji | (a normalization factor such as p · n can be neglected),
then a strategy that systematically beats random input assignment consists
in ranking the nodes according to the ratio rw (i) = wout (i)/win (i) and placing inputs on the nodes with highest rw . In Fig. 3(a) the λmin (Wm ) of this
driver node placement strategy is compared with a random selection. If for
full graphs the improvement is minimal, as the graphs become sparser it
increases, and for ER networks with p = 0.01 the λmin (Wm ) obtained by
controlling nodes with high rw is more than twice that of the random choice
of controls, see Fig. 3(b). As can be seen in Fig. S3, all measures of input

6

energy show a qualitatively similar improvement. When the topology of the
network renders the values of rw more extreme, like when direct scale-free
(SF) graphs are chosen, with indegree exponent bigger than outdegree exponent [5], see Fig. S4, then the improvement in choosing driver nodes with
high rw becomes much more substantial, even of orders of magnitude bigger
than a random selection, see Fig. 3(b) and Fig. S5 for more details.
What we deduce from such results is that once the technical issues associated with minimal controllability can be neglected, a general criterion for
controlling a network with a limited input cost is to drive the nodes having
the maximal disembalance between their weighted outdegree and indegree.
Notice that our computation of weighted out/indegrees considers the total
sum of weights in absolute value. When signs are taken into account in
computing win and wout , then no significant improvement over random input
placement is noticeable. This is connected to the quadratic nature of the
Gramian.
It is worth emphasizing that for a dynamical system the concept of driver
node is not intrinsic, but basis-dependent. In fact, just like the idea of
adjacency matrix of a graph is not invariant to a change of basis in state
space, so inputs associated to single nodes (i.e., to single state variables) in
the original basis become scattered to all variables in another representation
of the same system, see Fig. 4(a). If we take a special basis in which the modes
are decoupled (for instance the Jordan basis), then the contribution of the
nodes to the modes (i.e., the eigenvectors of A) provide useful information
for the investigation of minimum input energy problems. The topic is closely
related to the so-called participation factors analysis in power networks [30].
Also quantities like win and wout are basis-dependent and become nearly
equal for instance if in (1) we pass to a Jordan basis. On the contrary, the
eigenvalues of A are invariant to a change of basis. Hence as a general rule, the
control energy considerations that are consequence of the time constants of
the system (like the dependence on the real part of the eigenvalues illustrated
in Fig. 2) are “more intrinsic” than those that follow from the particular basis
representation we have available for a network.
Real part of the eigenvalues and controllability with bounded controls. From what we have seen above, the control energy is least when the
real part of the eigenvalues tends to vanish. In the limit case of all eigenvalues of A being purely imaginary, we recover a known result from control

7

theory affirming that controllability from any xo to any xf can be achieved
in finite time by means of control inputs of bounded amplitude. As a matter
of fact, an alternative approach used in control theory to take into account
the control cost of a state transfer is to impose that the amplitude of the
input stays bounded for all times (rather than the total energy), and to
seek for conditions that guarantees controllability with such bounded controls [6, 15, 18]. Assume u ∈ Ω, with Ω a compact set containing the origin,
for instance Ω = [−1, 1]m , where m is the number of control inputs. The
constraint u ∈ Ω guarantees that we are using at all times a control which
has an energy compatible with the physical constraints of our network. The
consequence is, however, that reaching any point in Rn may require a longer
time, or become unfeasible. In particular a necessary and sufficient condition
for any point to be reachable from 0 in finite time when u ∈ Ω is that no
eigenvalue of A has a negative real part, see SI. This is clearly connected with
our previous considerations on the reachability problem without bounds on
u: when all modes of A are unstable then the input energy required to reach
any state from 0 is low (Fig. 1(a)) and becomes negligible for sufficiently
long time horizons. On the contrary, transferring any state to 0 in finite
time with u ∈ Ω is possible if and only if no eigenvalue of A has a positive
real part. Also in this case the extra constraints on the input amplitude
reflects the qualitative reasoning stated above and shown in Fig. 1(b). Also
in the bounded control case, considering a generic transfer from any state
xo to any other state xf means combining the two scenarios just described:
formally a system is completely controllable from any xo to any xf in finite
time and with bounded control amplitude u ∈ Ω if and only if all eigenvalues
of A have zero real part, see SI for the details. The findings discussed above
for u unbounded are completely coherent with this alternative approach to
“practical controllability”.
Systems with purely imaginary eigenvalues: the case of coupled
harmonic oscillators. A special case of linear system with purely imaginary eigenvalues is a network of n coupled harmonic oscillators, represented
by a system of second order ODEs
M q̈ + Kq = Bu

(4)

where M = M T > 0 is the inertia matrix, K = K T > 0 is the stiffness
matrix, typically of the form K = Kd + L, with Kd > 0 diagonal and L a
8

Laplacian matrix representing the couplings. In (4) the controls are forces,
and the input matrix B contains elementary vectors in correspondence of the
controlled nodes. The state space representation of (4) is
ẋ = Ao x + Bo u

(5)

with


Mq
x=
∈ R2n ,
M q̇


Ao =


I
0
,
−KM −1 0


and

Bo =


0
.
B

The system (5) has purely imaginary eigenvalues equal to ±iωj , j = 1, . . . , n,
where ωj are
frequencies of the oscillators. If Ω2 = diag(ω12 , . . . , ωn2 )
 1the natural

n
and Ψ = ψ . . . ψ is the matrix of corresponding eigenvectors, then in the
so-called modal basis the oscillators are decoupled and one gets the state
space representation
ż = A1 z + B1 u




(6)
0
I
0
=
z+
u
2
T
−1
−Ω 0
Ψ M B

 −1
0
Ψ
x. See SI for the details. When a system has
where z =
0 Ψ−1
purely imaginary eigenvalues, the finite time Gramian diverges as tf → ∞.
However, in the modal basis (6) the Gramian is diagonally dominant and
linear in tf , hence as tf grows it can be approximated by a diagonal matrix
which can be computed explicitly [3]:
 (ψ1 )2

j




n
X βj 

Wz (tf ) ≈
2 
2Mj 
j=1




ω12

..






 tf ,





.
(ψjn )2
2
ωn

(ψj1 )2

..

.

(7)

(ψjn )2
where βj = 1 if the j-th input is present and 0 otherwise. Using (7), the
three measures of control energy adopted in this paper give rise to simple
strategies for minimum energy driver nodes placement, which in some cases
9

can be computed exactly for any n (for instance for the metric tr(Wz ), see
SI). Fig. 4 shows that such strategies are always beating a random driver
node placement, often by orders of magnitude.
Also wout /win is still a good heuristic for driver node placement strategy.
This can be understood by observing that the model (4) is symmetric hence
for it in- and out-degrees are identical. However, since Ao has rows rescaled
by M −1 , wout is affected
directly: when the inertia Mii is big, the correP
sponding wout (i) = nj=1 Kji /Mii is small and viceversa. No specific effect is
instead induced on win . In the representation (5), selecting nodes according
to wout /win means placing control inputs on the lighter masses, see Fig. 4 (d).
When the harmonic oscillators are decoupled (L = 0) then m < n means controllability is lost, but nevertheless the least control energy of the m inputs
is indeed obtained when driving the oscillators of least inertia. A weak (and
sparse) coupling allows to recover controllability, while the least inertia as optimal driver node strategy becomes suboptimal. When the coupling becomes
stronger (for instance when the coupling graph is more connected) then the
inertia of an oscillator is less significant as a criterion for selection of driver
nodes: the modes of the system are now spread throughout the network and
no longer localized on the individual nodes. As shown in Fig. S6, in a fully
connected network of harmonic oscillators, driver node strategies based on
wout /win and on tr(Wz ) tend to perform considerably worse for the other
measures (λmin (Wz ) and tr(Wz−1 )), while for a sparse graph (here ER graphs
with p = 0.05), of the three explicit optimal driver node placement strategies
available in this case, tr(Wz ) has a high overlap with wout /win , see Fig. 4
(b), while the other two tend to rank controls in somewhat different ways.
Given that in this case we have three strategies that are (near) optimal for
the chosen measure of control energy, the dissimilarity of the node rankings
of these three strategies means that the driver node placement problem is
heavily dependent on the way control energy is quantified.
Controlling vibrations of polyatomic molecules. Coupled harmonic
oscillators are used in several applicative contexts, for instance in the active
control of mechanical vibrations [12, 19] or in that of flexible multi-degree
structures (aircrafts, spacecrafts, see [17]) where our controllability results
can be applied straightforwardly (and compared with analogous/alternative
approaches described for instance in [3, 14, 17, 19, 32, 39, 42]). Another
context in which they are used is in controlling molecular vibrations of poly-

10

atomic molecules [43, 48]. The assumption of harmonicity is valid in the
regime of small oscillations near equilibrium, in which the potential energy
is approximated well by the parabola V (q) = 12 q T Kq (here q is the quantum
expectation value of the displacement from the equilibrium position and only
vibrational degrees of freedom are considered1 ). The node-driving setting described by (5) corresponds here to controlling the vibrations of single bonds,
as in mode-selective chemistry [9, 36, 48]. In the modal basis (6), the free
evolution of the modes (i.e., A1 ) is decoupled, but not the input matrix B1 ,
see Fig. 4(a). The B1 terms in (6) specify how the energy exciting a specific
vibrational bond propagates through the entire molecule and affects all its
modal frequencies. The methods developed above for quantifying the control
energy are applicable also to this context. In particular the Gramian can be
used to estimate the energy spreading of a monochromatic input field: for
the j-th bond it is proportional to the j-th row of Ψ (i.e., it consists of the
j-th components of all eigenvectors ψ 1 , . . . , ψ n ).
The basic principle we have adopted so far (minimize the input energy,
here the fluence of the pumping field) implicitly favours the selection of inputs
having a broad “coverage” in the modal basis, or, said otherwise, favours the
intramolecular spreading of vibrational energy to the entire molecule. This
is clearly not the main task of selective-mode laser chemistry, which on the
contrary aims at keeping as much energy as possible concentrated on specific
bonds or modes. Given that the power of the laser field is not a limiting
factor, a control problem that can be formulated using the tools developed in
this paper is for instance to choose m monochromatic excitations selectively
tuned on the stiffness constants of m bonds (i.e. for us certain columns
of Bo ) so as i) to guarantee controllability; ii) to maximize the power at
a certain mode ωk (representing for example a conformational change that
one wants to achieve in the molecule [9]). For the diagonal Gramian
Pm(7), this
amounts to choosing the indexes j1 , . . . , jm ∈ {1, . . . , n} such that `=1 (ψjk` )2
is maximized, a problem which can be solved exactly once Ψ is known. Notice
that even when energy spreads through the bonds because of the couplings,
it is in principle possible to “refocus” it towards a single bond using the
dynamics of the system. In the linear regime, a formula like (3) can be used
to compute explicitly the controls needed to refocus it on a specific bond,
1
For a molecule of n atoms the number of independent degrees of freedom is 3n − 5
i.e., 3 for each atom, minus the coordinates of the center of mass. All masses and stiffness
constants should be rescaled relative to it.

11

corresponding for instance to qf having a single non-zero component. This
does not require to solve an optimal control problem, as proposed for instance
in [34, 33].
Minimum energy control of power grids. In the linear regime, power
grids can be modeled as networks of weakly damped coupled harmonic oscillators [13]. The so-called swing equation corresponds in fact to the following
network of damped and coupled harmonic oscillators
M q̈ + Dq̇ + Kq = Bu,

(8)

where D is the matrix of dampings which we assume to be proportional, that
is, that in the modal basis D1 = ΨT M −1 DM −1 Ψ is diagonal. In the state
space representation (5), one gets then




0
I
0
x+
u,
(9)
ẋ =
−KM −1 −DM −1
B
while in the modal basis




0
I
0
ż =
z+
u.
−Ω2 −D1
ΨT M −1 B

(10)

For weak damping, the driver node selection strategies illustrated above can
be applied to the model (9) and so can the method based on wout /win . We
have investigated the minimum energy control of several power grids listed in
Table S1, varying the dampings across several orders of magnitude, see Fig. 5
(a). As expected, for all of them the energy required to achieve controllability
increases as the real part of the eigenvalues moves away from the imaginary
axis, see Fig. 5 (b) and Figs. S7-S10. All strategies still beat a random driver
node placement, even those based on the Gramian (7), formally valid only
for undamped dynamics.

12

SUPPLEMENTARY MATERIAL

4
4.1

Methods
Control energy: finite-time horizon formulation

Consider a linear system
ẋ = Ax + Bu

(S11)

where x ∈ Rn is the state vector, A ∈ Rn×n is the state update matrix,
B ∈ Rn×m is the input matrix, and u is the m-dimensional input vector. The
reachable set of (S11) in time tf from xo is the set
Rtf (xo ) = {x ∈ Rn s. .t. ∃ u : [0, tf ] → Ω s. t. φ(tf , u, xo ) = x}
where φ(t, u, xo ) is the solution of (S11) at time t with input u and Ω is the
admissible set of control inputs, here Ω = Rm .
The system (S11) is reachable (or controllable from the origin [2]) in time
tf if any xf ∈ Rn can be reached from 0 by some control u ∈ Ω in time tf ,
i.e. if Rtf (0) = Rn . It is said controllable to the origin if any xo ∈ Rn can
be brought to 0 by some control u ∈ Ω in time tf . The system (S11) is said
completely controllable in time tf if Rtf (xo ) = Rn for any xo ∈ Rn .
Finite-time Gramians. The time-tf reachability (or controllability from
0) Gramian is the symmetric matrix
Z tf
T
eAτ BB T eA τ dτ,
Wr (tf ) =
(S12)
0

while the time-tf controllability to 0 Gramian (normally called the controllability Gramian, [2]) is
Z tf
T
Wc (tf ) =
e−Aτ BB T e−A τ dτ.
(S13)
0

The two Gramians are positive definite whenever (A, B) is controllable, and
are related by
T
Wr (tf ) = eAtf Wc (tf )eA tf .
13

Similarly, for their inverses,
Wr−1 (tf ) = e−A

Tt
f

Wc−1 (tf )e−Atf .

Finite-time control energy for state transfer. The transfer of the state
from any xo into any other xf in time tf 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
Z
tf

ku(τ )k2 dτ.

E(tf ) =

(S14)

0

Such control can be computed explicitly [2] as
u(t) = B T eA

T (t −t)
f

Wr−1 (tf )(xf − E Atf xo ),

t ∈ [0, tf ]

(S15)

and the corresponding transfer cost as
E(tf ) = (xf − eAtf xo )T Wr−1 (tf )(xf − eAtf xo ).

(S16)

The various metrics that have been proposed in the recent and old literature to quantify the energy needed for state transfer between any two states
xo and xf are in fact all based on the Gramian [24]:
1. λmin (Wr ) = λmax (Wr−1 ): the min eigenvalue of the Gramian (equal
to the max eigenvalue of Wr−1 ) is a worst-case metric, estimating the
energy required to move along the direction which is most difficult to
control.
2. tr(Wr ): the trace of the Gramian is inversely proportional to the average energy required to control a system.
3. tr(Wr−1 ): the trace of the inverse of the Gramian is proportional to the
average energy needed to control the system.
Minimizing the control energy means maximizing the first and second
measure or minimizing the third. To be more precise on how the control energy is formed, we have to split the state transfer energy (S16) into subtasks:
1. xo = 0 (reachability problem)
=⇒

Er (tf ) = xTf Wr−1 (tf )xf ;
14

2. xf = 0 (controllability to 0 problem)
=⇒

Ec (tf ) = xTo eA

Tt
f

Wr−1 (tf )eAtf xo = xTo Wc−1 (tf )xo .

In particular, both Wr (tf ) and Wc (tf ) enter into the cost function. In particular, to quantify the amount of control energy of these problems we need
to compute the inverse of Wr (tf ) and Wc (tf ).
Let us look at how the stability/instability of the eigenvalues influences
the two costs Er (tf ) and Ec (tf ).
• 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 xo → xf requires to take into account the “difficult” directions of
both cases.

4.2

Control energy: infinite-time horizon formulation

When tf → ∞, then E(tf ) converges (or diverges) to a quantity
Z ∞
E=
ku(τ )k2 dτ,

(S17)

0

and so do Er (tf ) and Ec (tf ).
When tf → ∞, both Gramians become infinite-time integrals, which may
be convergent or divergent, depending on the modes of A. If A stable, then
Z ∞
T
Wr =
eAτ BB T eA τ dτ
(S18)
0

15

exists finite and it is positive definite if (A, B) controllable. If instead A is
antistable, then it is
Z ∞
T
e−Aτ BB T e−A τ dτ
(S19)
Wc =
0

to exist finite and positive definite when (A, B) controllable. In the mixed
eigenvalues cases the two expression (S18) and (S19) both diverge.
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, limt→∞ x(t) = 0 for all
xo , meaning that the controllability to 0 problem can be solved with
zero energy Ec = 0. Furthermore, the integral (S12) converges to (S18),
whose value can be computed solving the following Lyapunov equation:
AWr + Wr AT + BB T = 0.

(S20)

Such a solution always exists and it is Wr > 0 (positive definite) if the
pair (A, B) is controllable.
• When instead A is antistable, then the integral (S12) diverges as tf →
∞. Hence the solution with u = 0 is no longer feasible as all modes
are unstable (and diverge as soon as xo 6= 0), meaning that to find a
minimizer of (S17) 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)
P (−A) + (−AT )P + P BB T P = 0.

(S21)

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 exist)
(−A)L + L(−AT ) + BB T = 0,
16

(S22)

and then setting P = L−1 . It can be verified directly that the controllability Gramian Wc in (S19) is one such solution L. Correspondingly
we obtain P = Wc−1 . From the theory of linear-quadratic regulators
(in particular [40], Ch. 10) the feedback controller
u = −B T P x(t)

(S23)

guarantees stability of the closed-loop system
ẋ = (A − BB T P )x
i.e., A − BB T P is a stable matrix. The feedback law (S23) also minimizes the input energy (S17) which is equal to Ec = xTo Wc−1 xo .
• 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 [49]. Split B and x(t) accordingly:
 


 
x1
A1 0
B1
Re[λ(A1 )] < 0
x=
, A=
, B=
,
x2
0 A2
B2
Re[λ(A2 )] > 0
(S24)
In the infinite time horizon, the u = 0 control steers optimally the
x1 subvector, while for the x2 part a feedback controller provides the
energy-minimizing solution. From (S21) we obtain that the ARE has
solution


0 0
P =
0 P2
where P2 solves the ARE for the (A2 , B2 ) subsystem. Hence the control
input


0
T
T 0
u = −BB P = −BB
−1
0 W2,c
achieves a transfer to the origin with minimal energy cost equal to
−1
Ec = xTo P xo = xT2,o W2,c
x2,o . Furthermore, combining (S20) and (S22),
we have that for (S24) the following two Lyapunov equations must hold
simultaneously:
A1 W1,r + W1,r AT1 + B1 B1T = 0
(−A2 )W2,c + W2,c (−AT2 ) + B2 B2T = 0
17

(S25a)
(S25b)

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 + AT P + P BB T P = 0

(S26)

or, equivalently, P = K −1 with K solution of the Lyapunov equation
AK + KAT + BB T = 0.

(S27)

A is stable, hence K > 0 solving (S27) and P > 0 solving (S26) always
exist. The resulting closed loop matrix A − BB T P must be antistable.
From (S20) and (S27) it can also be K = Wr , 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 xf when initialized in xo = 0; this amount of energy
is negligible in the infinite time horizon considered here). Since −A
stable, a Lyapunov equation like (S22) holds, with solution L = Wc .
• When A has eigenvalues with both positive and negative real part and
no purely imaginary eigenvalues, then a decomposition like (S24) can
be obtained through a change of basis. The complete ARE has now
solution

  −1 
W1,r 0
P1 0
P =
=
,
0 0
0
0
and the controller achieving the transfer with minimal energy is
 −1 
0
T W1,r
,
u = −BB
0
0
for an amount of energy equal to
−1
Er = xT1,f W1,r
x1,f

The decomposition (S24) also in this case induce a pair of Lyapunov
equations identical to (S25).
18

4.3

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, Wm and which
gathers the directions difficult to control of both the reachability and the
controllability to 0 problems.
Infinite-time horizon mixed Gramian Assume that the spectrum of
A contains k, 0 6 k 6 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 (S24):


Ā1 0
= V AV −1
(S28)
0 Ā2
and, correspondingly,



B̄1
=VB
B̄2

with Re[λ(Ā1 )] < 0 and Re[λ(Ā2 )] > 0. In the new basis, the two Lyapunov
equations (S25) hold, which can be rewritten as


with

T

T 

0
B̄1 B̄1T
Ā1
0
Ā1
0
=0
+
W̄ + W̄
0
B̄2 B̄2
0 −Ā2
0 −Ā2

(S29)



W̄1,r
0
W̄m =
0
W̄2,c

the mixed Gramian. Following [49], the expression of the mixed Gramian in
the original basis is Wm = V −1 W̄m V −T . By construction, the mixed Gramian
matrix Wm always exists and summarizes the infinite-horizon contribution
of the stable eigenvalues to the reachability problem and of the unstable
eigenvalues to the controllability to 0 problem.
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

19

diagonal blocks, (S28), this is given by


W̄1,r (tf )
0
W̄m (tf ) =
0
W̄2,c (tf )
where W̄1,r (tf ) and W̄2,c (tf ) are the equivalent of (S12) and (S13) for the two
subsystems (Ā1 , B̄1 ) and (Ā2 , B̄2 ). An equation like (S29) has no coupling
terms between the two subsystems (i.e., terms of the form B̄1 B̄2 ). These
terms disappear asymptotically, but transiently they give a contribution,
hence the finite-time formulation of W̄m (tf ) is only an approximation. In the
original basis, Wm (tf ) = V −1 W̄m (tf )V −T , and the input energy is



 T
−1 x1,f
T
.
Em (tf ) = x1,f x2,o Wm (tf )
x2,o
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, or shifting the eigenvalues of
A by adding a diagonal term (see Fig. 2). Equivalently, all eigenvalues can
be made unstable, and the controllability to 0 problem considered.

4.4

Controllability with bounded controls

Consider the system (S11). Assume u ∈ Ω, where Ω is a compact set of Rm
containing the origin in its interior. Assume (A, B) is controllable. Then we
have the following, see [18, 6] and [35], p. 122.
• A necessary and sufficient condition for the origin to be steered to any
point of Rn 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 Rn 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 xo to any point xf ) in finite time is that all eigenvalues have
zero real part.
20

4.5

Control of coupled harmonic oscillators

A network of n coupled harmonic oscillators can be written as a system of
second order differential equations
Mi q̈i + (ki +

n
X
j=1

kij )qi −

n
X

kij qj = βi ui ,

i = 1, . . . , n,

(S30)

j=1

where Mi > 0 is the mass of the i-th oscillator, ki > 0 its stiffness, kij > 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, (S30) can be rewritten as
M q̈ + Kq = Bu,

(S31)

where M = M T = diag(Mi ) > 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 (S31) have
the form q = φeiωt in correspondence of the pairs ωj and φj , j = 1, . . . , n,
that are the solutions of the generalized eigenvalues/eigenvector equation
(−ω 2 M + K)φ = 0.

(S32)

The ωj are called the natural frequencies of (S31). Denote


Φ = φ1 . . . φn
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 q1 = Φ−1 q, M1 = ΦT M Φ and
K1 = ΦT KΦ are both diagonal matrices, hence
M1 q̈1 + K1 q1 = ΦT Bu
has decoupled dynamics (but coupled inputs).
The state space representation of (S31) is 2n dimensional. If

 
M 0
q
x=
,
0 M q̇
then


ẋ = Ao x + Bo u =




I
0
0
x+
u.
−KM −1 0
B
21

(S33)

(S34)

In terms of the state space model (S34), the eigenvalues are λj = ±iωj ,
j = 1, . . . , n, of eigenvectors


ψj
vj =
ωj ψ j
where ψ j = M φj , from which the purely oscillatory nature of Ao is evident.
If we denote


ω12


ω22


2
Ω =
,
.
.


.
2
ωn
then, from (S32), M −1 KΦ = ΦΩ which implies
Ω2 = Φ−1 M −1 KΦ
= Φ−1 M −1 Φ−T ΦT KΦ
= M1−1 K1 .
If Ψ = M Φ, the state space representation in the modal basis
 −1

Ψ
0
z = Tx =
x
0 Ψ−1
is given by

ż = A1 z + B1 u =




0
0
I
z+
u.
−Ω2 0
Ψ−1 B

(S35)

From M1 = ΨT M −1 Ψ, one gets Ψ−1 B = M1−1 ΨT M −1 B.
A key advantage of the modal representation is that the Gramian of the
pair (A1 , B1 ) can be computed explicitly. As a matter of fact, when the
eigenvalues are on the imaginary axis, the integral (S12) (or (S13)) diverge,
and hence the infinite-time Gramian cannot be computed. However, in the
modal basis z, Wz (tf ) (or more precisely Wz,r (tf )) is diagonally dominant,
and for tf sufficiently long it can be approximated by its diagonal terms.
These terms are computed explicitly in [3]:
 −1 T −1 T −1 −T
 (M1 Ψ M BB M ΨM1 )jj tf
for 1 6 j 6 n
2ωj
(Wz (tf ))jj =
−1 T
−1 BB T M −1 ΨM −T
M
Ψ
M
t
( 1
1 )jj f
for n + 1 6 j 6 2n.
2
22

If we assume that the mass matrix M is diagonal, then it is always possible
to choose Ψ so that M1 = I and K1 diagonal, by suitably rescaling the
eigenvectors ψ j . In this case






0
0
0
T
T
−1
0 B M Ψ =
,
B1 B1 =
ΨT M −1 B
0 ΨT M −1 BB T M −1 Ψ
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:
 
1
n
X βj  ψ.j  

1
n
.
ψ
.
.
.
ψ
ΨT M −1 BB T M −1 Ψ =


.
j
j
Mj2
j=1
ψjn
(only m of the n factors βj ∈ {0, 1} are nonzero) and its diagonal entries are


1 2
(ψ
)
n
j
X βj 

..
(S36)
diag(ΨT M −1 BB T M −1 Ψ) =
.
.
2 
Mj
n 2
j=1
(ψj )
Hence the expression for the Gramian in the modal basis is
P
1 2






Wz (tf ) ≈ 







βj (ψj )
n
j=1 Mj2 2ω12

..



.
n 2
βj (ψj )
2
j=1 Mj2 2ωn

Pn

1 2
βj (ψj )
j=1 Mj2 2

Pn

..

.
βj
j=1 Mj2

Pn

Notice the linearity in tf , meaning that all components diverge to ∞ with
the same speed when tf → ∞. 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

23







 tf .






(ψ n )2
j

2

inputs according to λmin (Wz ) amounts to solving the following MILP maxmin problem:
max min
βj

i

n
X
βj i 2
(ψj )
2
M
j
j=1

subject to

n
X

βi = m

j=1

βj ∈ {0, 1}
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 tr(Wz ).
Input:
yi =

1 h (ψi1 )2
...
2
Mi2 2ω1

(ψin )2
2
2ωn

(ψi1 )2
2

...

(ψin )2
2

i

,

i = 1, . . . , n

1. Choose î = argmaxi (−ky i k∞ ) = argmaxi (mink yki )
• y s = y î
• I = {1, 2, . . . , n} \ { î }
• O = { î }
2. For c = 2, 3, . . . , m
• compute ĵ = argmaxj∈I (−ky s + y j k∞ )
• y s = y s + y ĵ
• I = I \ { ĵ }
• O = O ∪ { ĵ }
Output: O

24

If instead we choose to maximize the tr(Wz ), then we get
n X
n
X
βj (ψji )2
max
βj
Mj2 2(1 + ωi2 )
i=1 j=1

subject to

n
X

βj = m

j=1

βj ∈ {0, 1}
Since tr(Wz ) is linear in the βj , this is a linear optimization problem, hence
solvable exactly and efficiently for any n. Finally, also for the minimization
of tr(Wz−1 )
min tr(Wz−1 )
βj

subject to

n
X

βj = m

j=1

βj ∈ {0, 1}
an efficient heuristic can be set, as outlined in Algorithm 2.

25

Algorithm 2 Driver node placement that minimizes tr(Wz−1 ).
Input:
yi =

1 h (ψi1 )2
...
2
Mi2 2(1+ω1 )

1. Choose î = argmaxi

(ψin )2
2)
2(1+ωn

i

,

i = 1, . . . , n

Pn

1
k=1 yki

• y s = y î
• I = {1, 2, . . . , n} \ { î }
• O = { î }
2. For c = 2, 3, . . . , m
• compute ĵ = argmaxj∈I

Pn

1
k=1 y s +y j
k
k

• y s = y s + y ĵ
• I = I \ { ĵ }
• O = O ∪ { ĵ }
Output: O

Looking at an expression like (S36), it is possible to understand what kind
of behavior yields good controllability properties to certain driver nodes. For
instance when measuring according to λmin (Wz ), from (S36), 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 ψj1 , . . . , ψjn are nonvanishing and possibly all as large as
possible. Since Wz 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 {ψik }k=1...,n and {ψjk }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 [12] for the details. In that case
26

ΨT = Ψ−1 i.e., the “coverage” discussed here is given by the left eigenvectors
of (S32), condition sometimes considered in the literature [14].
Another basis
  for the state space that can be used in place of (S33) is
q
given by x̃ =
. With this choice, the state space realization is
q̇




I
0
0
˙x̃ = Ão x̃ + B̃o u =
x̃ +
u.
(S37)
−M −1 K 0
M −1 B
It is straightforward to verify that in this basis the roles of win and wout
are exchanged, hence a criterion for driver node selection becomes ranking
according to win /wout (instead of wout /win ).

5

Datasets

The power networks used in the paper are listed in Table S1. All nodes are
treated equally, regardless of their function as generators or loads in the real
grid.

Network type

nodes

edges

source

North EU power grid
IEEE 300 test grid
French power grid
USA power grid

236
300
1888
4941

320
409
2308
6591

[23]
https://www.ee.washington.edu/research/pstca/
[16]
[44]

Table S1: Power grids used in this study.

References
[1] Stefano Allesina and Si Tang. The stability–complexity relationship at
age 40: a random matrix perspective. Population Ecology, 57(1):63–75,
2015.
[2] P.J. Antsaklis and A.N. Michel. Linear Systems. Birkhäuser Boston,
2005.
[3] Ami Arbel. Controllability measures and actuator placement in oscillatory systems. International Journal of Control, 33(3):565–574, 1981.
27

[4] N. Bof, G. Baggio, and S. Zampieri. On the Role of Network Centrality
in the Controllability of Complex Networks. ArXiv e-prints, September
2015.
[5] Béla Bollobás, Christian Borgs, Jennifer Chayes, and Oliver Riordan.
Directed scale-free graphs. In Proceedings of the Fourteenth Annual
ACM-SIAM Symposium on Discrete Algorithms, SODA ’03, pages 132–
139, Philadelphia, PA, USA, 2003. Society for Industrial and Applied
Mathematics.
[6] R. F. Brammer. Controllability in linear autonomous systems with positive controllers. SIAM J of Control, 10:339–353, 1972.
[7] Yu-Zhong Chen, Le-Zhi Wang, Wen-Xu Wang, and Ying-Cheng Lai.
Energy scaling and reduction in controlling complex networks. Royal
Society Open Science, 3(4), 2016.
[8] Sean P. Cornelius, William L. Kath, and Adilson E. Motter. Realistic
control of network dynamics. Nat Commun, 4, 06 2013.
[9] Brian C. Dian, Asier Longarte, and Timothy S. Zwier. Conformational
dynamics in a dipeptide after single-mode vibrational excitation. Science, 296(5577):2369–2373, 2002.
[10] Jin Ding, Yong-Zai Lu, and Jian Chu. Studies on controllability of
directed networks with extremal optimization. Physica A: Statistical
Mechanics and its Applications, 392(24):6603 – 6615, 2013.
[11] Jianxi Gao, Yang-Yu Liu, Raissa M D’Souza, and Albert-László
Barabási. Target control of complex networks. Nature communications,
5, 2014.
[12] Wodek K. Gawronski. Dynamics and Control of Structures. A Modal
Approach. Mechanical Engineering Series. Springer-Verlag, New York,
1998.
[13] L.L. Grigsby. Power System Stability and Control. The Electric Power
Engineering Hbk, Second Edition. CRC Press, 2007.
[14] A. M.A. Hamdan and A. H. Nayfeh. Measures of modal controllability
and observability for first- and second-order linear systems. Journal of
Guidance, Control, and Dynamics, 12:421–428, 1989.
28

[15] D.H. Jacobson. Extensions of Linear-Quadratic Control, Optimization
and Matrix Theory, volume 133 of Mathematics in Science and Engineering. Academic Press, London, 1977.
[16] C. Josz, S. Fliscounakis, J. Maeght, and P. Panciatici. AC Power Flow
Data in MATPOWER and QCQP Format: iTesla, RTE Snapshots, and
PEGASE. ArXiv e-prints, March 2016.
[17] J.L. Junkins and Y. Kim. Introduction to Dynamics and Control of
Flexible Structures. AIAA education series. American Institute of Aeronautics & Astronautics, 1993.
[18] E.B. Lee and L. Markus. Foundations of Optimal Control Theory. R.E.
Krieger Publishing Company, 1986.
[19] S. Leleu, H. Abou-Kandil, and Y. Bonnassieux. Piezoelectric actuators and sensors location for active control of flexible structures. IEEE
Transactions on Instrumentation and Measurement, 50(6):1577–1582,
Dec 2001.
[20] Guoqi Li, Wuhua Hu, Gaoxi Xiao, Lei Deng, Pei Tang, Jing Pei, and
Luping Shi. Minimum-cost control of complex networks. New Journal
of Physics, 18(1):013012, 2016.
[21] Y.-Y. Liu and A.-L. Barabási. Control Principles of Complex Networks.
ArXiv e-prints, August 2015.
[22] Yang-Yu Liu, Jean-Jacques Slotine, and Albert-Laszlo Barabasi. Controllability of complex networks. Nature, 473(7346):167–173, 2011.
[23] Peter J. Menck, Jobst Heitzig, Jürgen Kurths, and Hans
Joachim Schellnhuber. How dead ends undermine power grid stability. Nat Commun, 5, 06 2014.
[24] P.C. Müller and H.I. Weber. Analysis and optimization of certain qualities of controllability and observability for linear dynamical systems.
Automatica, 8(3):237 – 246, 1972.
[25] Jose C. Nacher and Tatsuya Akutsu. Analysis of critical and redundant
nodes in controlling directed and undirected complex networks using
dominating sets. Journal of Complex Networks, 2(4):394–412, 2014.
29

[26] Tams Nepusz and Tams Vicsek. Controlling edge dynamics in complex
networks. Nature Physics, 05 2012.
[27] A. Olshevsky. Minimal controllability problems. Control of Network
Systems, IEEE Transactions on, 1(3):249–258, Sept 2014.
[28] A. Olshevsky. Eigenvalue Clustering, Control Energy, and Logarithmic
Capacity. ArXiv e-prints, October 2015.
[29] F. Pasqualetti, S. Zampieri, and F. Bullo. Controllability metrics, limitations and algorithms for complex networks. IEEE Transactions on
Control of Network Systems, 1(1):40–52, March 2014.
[30] I. J. Perez-arriaga, G. C. Verghese, and F. C. Schweppe. Selective modal
analysis with applications to electric power systems, part i: Heuristic
introduction. IEEE Transactions on Power Apparatus and Systems,
PAS-101(9):3117–3125, Sept 1982.
[31] Justin Ruths and Derek Ruths. Control profiles of complex networks.
Science, 343(6177):1373–1376, 2014.
[32] Hamid Reza Shaker and Maryamsadat Tahavori. Optimal sensor and actuator location for unstable systems. Journal of Vibration and Control,
2012.
[33] Shenghua Shi and Herschel Rabitz. Optimal control of selective vibrational excitation of harmonic molecules: Analytic solution and restricted forms for the optimal fields. The Journal of Chemical Physics,
92(5):2927–2937, 1990.
[34] Shenghua Shi, Andrea Woody, and Herschel Rabitz. Optimal control of
selective vibrational excitation in harmonic linear chain molecules. The
Journal of Chemical Physics, 88(11):6870–6883, 1988.
[35] Eduardo D. Sontag. Mathematical Control Theory: Deterministic Finite
Dimensional Systems (2Nd Ed.). Springer-Verlag New York, Inc., New
York, NY, USA, 1998.
[36] Peter F. Staanum, Klaus Hojbjerre, Peter S. Skyt, Anders K. Hansen,
and Michael Drewsen. Rotational laser cooling of vibrationally and
translationally cold molecular ions. Nat Phys, 6(4):271–274, 04 2010.
30

[37] T. H. Summers, F. L. Cortesi, and J. Lygeros. On submodularity and
controllability in complex dynamical networks. IEEE Transactions on
Control of Network Systems, 3(1):91–101, March 2016.
[38] Jie Sun and Adilson E. Motter. Controllability transition and nonlocality
in network control. Phys. Rev. Lett., 110:208701, May 2013.
[39] M. Tarokh. Measures for controllability, observability and fixed modes.
IEEE Transactions on Automatic Control, 37(8):1268–1273, Aug 1992.
[40] H. Trentelman, A.A. Stoorvogel, and M. Hautus. Control Theory for
Linear Systems. Communications and Control Engineering. Springer
London, 2012.
[41] V. Tzoumas, M. A. Rahimian, G. J. Pappas, and A. Jadbabaie. Minimal
actuator placement with bounds on control effort. IEEE Transactions
on Control of Network Systems, 3(1):67–78, March 2016.
[42] Marc van de Wal and Bram de Jager. A review of methods for input/output selection. Automatica, 37(4):487 – 510, 2001.
[43] W. S. Warren, H. Rabitz, and M. Dahleh. Coherent control of quantum
dynamics: The dream is alive. Science, 259:1581, 1993.
[44] Duncan J. Watts and Steven H. Strogatz. Collective dynamics of smallworld networks. Nature, 393(6684):440–442, 06 1998.
[45] Gang Yan, Jie Ren, Ying-Cheng Lai, Choy-Heng Lai, and Baowen Li.
Controlling complex networks: How much energy is needed? Phys. Rev.
Lett., 108:218703, May 2012.
[46] Gang Yan, Georgios Tsekenis, Baruch Barzel, Jean-Jacques Slotine,
Yang-Yu Liu, and Albert-Laszlo Barabasi. Spectrum of controlling and
observing complex networks. Nat Phys, 11(9):779–786, 09 2015.
[47] Zhengzhong Yuan, Chen Zhao, Zengru Di, Wen-Xu Wang, and YingCheng Lai. Exact controllability of complex networks. Nat Commun, 4,
09 2013.
[48] Ahmed H. Zewail. Laser selective chemistry: is it possible? Physics
Today, 33(11):27–33, 1980.
31

[49] Kemin Zhou, Gregory Salomon, and Eva Wu. Balanced realization and
model reduction for unstable systems. International Journal of Robust
and Nonlinear Control, 9(3):183–198, 1999.

32

Energy

A mixed

A stable

2

2

0
1

0
1
0
-1

x2

A antistable
2
0
1
0

0
0

-1

1

x1

-1

x2

-1

0

1

x1

-1

x2

-1

0

1

x1

(a)

Energy

A mixed

A stable

A antistable

2

2

2

0
1

0
1

0
1

0
-1

x2

0

0
0

-1

1

x1

-1

x2

-1

0

1

x1

-1

x2

-1

0

1

x1

(b)

Energy

| Re[λ (A)] | <2

| Re[λ (A)] | <0.5

5

0
1

0
1

0
1
0

x2

| Re[λ (A)] | <0.1
5

5

-1

0

0
0

-1

x1

1

x2

-1

-1

0

x1

1

x2

-1

-1

0

1

x1

(c)

Figure 1: Reachability and Controllability-to-0 problems. (a): The reachability
(or controllability-from-0) problem is difficult along the stable eigendirections of A
(red curves in the leftmost panel) and easy along the unstable ones (blue). This is
reflected in the surfaces of Er (tf ) = xTf Wr−1 (tf )xf shown in the 3 rightmost panels.
In particular, the reachability problem requires limited control energy when A is
antistable (rightmost panel). (b): The controllability-to-0 problem is difficult along
the unstable eigendirections of A (red) and easy along the stable ones (blue). The
input energy surfaces, Ec (tf ) = xTo Wc−1 (tf )xo , reflect these properties. The case of
A stable requires the least control energy. (c): The problem studied in this paper
is a mixture of the two cases, collecting the worst-case of both. When the real
part of the eigenvalues of A is squeezed towards the imaginary axis as in the lower
right panels of Fig. 2(a), the input energy reduces accordingly.
33

(a)

tr(Wm )/n

λ min (Wm )

0.6
0.4

10

tr(W-1
)/n
m

10 4

1
0.8

2

10 10

10 5

0.2
10 0

10 0

200

400

n. inputs

600

800

200

400

n. inputs

600

800

200

400

600

800

n. inputs

(b)

Figure 2: (a): Circular law and eigenvalue location. For a random matrix, the
circular law allows to obtain state matrices A with eigenvalues in prescribed locations, for instance in the unit disk (blue) or in the shifted unit disk (red and yellow)
by altering the diagonal of A. The elliptic law allows to squeeze the eigenvalue
location along one of the two axes of the complex plane (violet, green and cyan).
(b): Control energy for various metrics when the number of (randomly chosen)
inputs grows. The data show a mean over 100 realizations of dimension n = 1000
(for each realization 100 different edge weights assignments are considered). The
color code is as in (a). For all three metrics used to measure the control energy
−1 ) which should
(λmin (Wm ), tr(Wm ) which should both be maximized, and tr(Wm
be minimized), the performances are strictly a function of the position of the eigenvalues of A. The minimum of the control energy is achieved when the eigenvalues
have very small real part (cyan) and worsen with growing real part, following the
order: cyan, green, blue, violet.

34

SF, γ =3.14, γ

ratio of λ min (Wm )

0.1

λmin (Wm )

ER, p = 0.05
0.08
0.06
0.04

ratio w out /w in

0.02

out

=2.87

ratio of λ min (Wm )

in

102

2

ER, p = 0.01

1.5

ER, p = 0.05

101

full
1
0

200

400

600

n. inputs

random
100

200

300

400

500

600

100

n. inputs

0

100

200

300

400

500

600

n. inputs

(a)

(b)

(c)

Figure 3: Driver node placement strategy: ranking according to rw = wout /win .
(a): When the driver nodes are chosen according to the ratio rw , then all the
control energy measures improve with respect to a random node selection. Here
λmin (Wm ) is shown, the other energy measures are in Fig. S3. Measures are means
(and st. dev.) over 100 realizations of size n = 1000; for each realization 100
edge weight assignments are tested. (b): For ER networks, the improvement in
λmin (Wm ) increases with the sparsity of the graph (inset: zoomed comparison in
linear scale). For other topologies, like SF directed graphs with indegree exponent
γin = 3.14 and outdegree exponent γout = 2.87, the improvement is remarkably
more significant (two orders of magnitude, violet curve, see also Fig. S5 for more
details). (c): the ratio rw of ranked nodes is shown. For SF networks, the fraction
of nodes having high rw is much bigger than that of ER networks, and this leads
to the much better performances in terms of control energy. The shaded areas
represent the values of m tested in our computations.

35

(a)
10 0

10 15

0.02
tr(Wz )
tr(W-1
)
z
w out /w in

0.01

random

tr(W-1
)/n
z

0.015

tr(Wz)/n

λ min (Wz)

λ min (Wz )

10 -1

0.005

100

200

300

400

500

100

600

200

300

400

500

10 10

10 5

600

100

200

n. inputs

n. inputs

300

400

500

600

8

10

n. inputs

(b)

0.5

0

overlap with w in/w out

0.5

0
0

200

400

n. inputs

600

1

overlap with λ min (Wz )

0.5

0
0

200

400

n. inputs

(c)

600

2

w out (i)/win (i)

1

%

overlap with w out /w in

%

%

1

0

200

400

600

1.5
1
0.5
0
0

2

4

6

Mi

n. inputs

(d)

Figure 4: Driver node placement strategies for a network of coupled harmonic
oscillators. (a): The concept of driver node is basis dependent: when the basis
changes in state space (for instance we pass from (5) to (6)), the control inputs no
longer target a single node, but become spread across the entire state space (now
decoupled into non-interacting modes). (b): Comparison of different driver node
placement strategies for n = 1000 coupled harmonic oscillators. Shown are means
over 100 realizations (with 100 edge weights samples taken for each realization).
Red: driver node placement based on λmin (Wz ). Violet: placement based on
tr(Wz ). Green: placement based on tr(Wz−1 ). Cyan: placement based on wout /win .
Blue: random input assignment. All driver node placement strategies always beat
a random assignment, often by orders of magnitude. The green and red curves
give similar performances and so do the cyan and violet. Notice that for tr(Wz )
the violet curve gives the exact optimum. (c): Overlap in the node ranking of the
different driver node placement strategies. Color code is the same as in (b). The
only highly significant overlap is between wout /win and tr(Wz ), while λmin (Wz ) and
tr(Wz−1 ) correspond to different node ranking patterns. Notice that none of the
strategies orders nodes according to win /wout (mid panel). (d) Inverse correlation
between Mi and wout /win (correlation coefficient
around −0.75 in average).
36

(a)
×10 -3
12

10

λ min (W) placement

w out /w in placement

×10 -6
14

-2

0.05

10

-2

12

8

10

6

10 -15

4

50

100

0.04
0.03

λ min (W)

-10

λ min (W)

λ min (W)

10

10 -10
10 -14

0.02

50

100

0.01

2
40

60

n. inputs

80

100

10

10 -10

8
6

random placement

10 -5

10

-15

50

4

100

2
40

60

n. inputs

80

100

40

60

80

100

n. inputs

(b)

Figure 5: Minimum energy control of power grids for varying damping coefficients.
(a): The eigenvalues of the state space system (9) for the North EU power grid
[23] with uniformly distributed masses (hMi i = 10) and damping coefficients that
vary across 4 orders of magnitude. (b): Control energy for the metric λmin (Wr )
when the driver nodes are placed according to λmin (Wr ) (left panel), wout /win
(mid panel), or randomly (right panel). The values of λmin (Wr ) corresponding
to the 4 choices of damping made in (a) are shown in solid lines (same color
code as in (a)), while in dotted lines the values of λmin (Wz ) are shown (suitably
normalized to eliminate the explicit dependence from tf , see (7)). The insets show
the same quantities in log scale. Values are all averages over 100 realizations.
For all three driver node placement strategies, the performances worsen as the
damping is increased. Comparing the three panels, wout /win performs similarly to
λmin (Wr ), and both outperform a random placement by orders of magnitude.

37

(a)

3

2
1.5
1

10 4

tr(W-1
)/n
m

tr(Wm )/n

λ min (Wm )

2.5

10 2

0.5
400

n. inputs

600

800

10 5

10 0

10 0

200

10 10

200

400

n. inputs

600

800

200

400

600

800

n. inputs

(b)

Figure S1: Analogous of Fig. 2, but for networks with an ER topology with
edge probability p = 0.05. (a): Six different eigenvalues locations in the complex
plane for ER graphs of size n = 1000 and random edges weights. The circular
law and the elliptic law are still valid. (b): Control energy for various metrics
when the number of (randomly chosen) inputs grows. The data show a mean over
100 realizations of dimension n = 1000 (for each realization 100 different edge
weights assignments are considered). The color code is as in (a). For all three
metrics used to measure the control energy (λmin (Wm ), tr(Wm ) which should both
−1 ) which should be minimized), the performances are
be maximized, and tr(Wm
strictly a function of the position of the eigenvalues of A. The minimum of the
control energy is achieved when the eigenvalues have very small real part (cyan)
and worsen with growing real part, following the order: cyan, green, blue, violet.

38

5000
4000

0.25

tr(Wm )/n

λ min (Wm )

0.3

0.2
0.15

3000
2000

0.1
1000

0.05
20

100

200

20

∞

100

time

200

∞

200

∞

time
×10 7
2.5

λ max (Wm )/λ min (Wm )

120

tr(W-1
)/n
m

100
80
60
40
20
20

100

200

∞

time

2
1.5
1
0.5
20

100

time

Figure S2: Computing control energies in finite time and infinite time. For the
various measures of control energy considered in the paper (λmin (Wm ), tr(Wm )
−1 )), the plots show the profile in time when computations are performed
and tr(Wm
using Wm (tf ), for various values of tf . The value for tf = ∞ is also shown for
comparison. For this specific example (full network of size n = 1000 and m = 400
−1 )
controls), some measures converge much faster than others. For instance tr(Wm
achieves its infinite-time value extremely quickly, while tr(Wm ) converges very
slow. Also the condition number of Wm (i.e., λmax (Wm )/λmin (Wm ), lower right
panel) converges to its asymptotic value.

39

0.1
1500

tr(Wm )/n

λmin (Wm )

0.08
0.06
0.04

500

0.02

200

300

400

500

600

λmax(Wm )/λmin (Wm )

100

tr(W-1
)/n
m

10 6

10

10

1000

4

2

100

200

300

400

500

100

200

100

200

300

400

500

600

300

400

500

600

10 12

10 10

10

8

600

n. inputs

n. inputs

ratio of tr(W m )

ratio of λmin (Wm )

(a)

1.3

1.2

1.1

1

1.2

1.1

1
200

400

600

1
0.95
0.9
0.85
0.8
0.75
0

200

400

600

ratio of λmax(Wm )λmin (Wm )

0

ratio of tr(W -1
)
m

1.3

n. inputs

0

200

400

600

0

200

400

600

1
0.95
0.9
0.85
0.8
0.75

n. inputs

(b)

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 rw = wout /win (red) and the same measure
for random driver node assignments (blue). 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 Wm (i.e., λmax (Wm )/λmin (Wm ), lower right panel) improves.

40

cum. prob.

10 0

cum. prob.

10 0

ER, p=0.05
ER, p=0.01
SF, γ in =3.14, γ out =2.87

10 -5
10 0

10 1

10 -5
10 0

10 2

10 1

indegree

10 2

outdegree

(a)

3

full
ER, p=0.05
ER, p=0.01
SF

d out /d in

2.5
2
1.5
1
10

1

10

2

10

3

nodes

(b)

Figure S4: Degree distributions of the various networks used in this study. (a):
Average indegree and outdegree of the ER networks (with p = 0.05 and p = 0.01)
and SF networks. The SF networks are generated using the algorithm of [5],
with indegree exponent γin = 3.14 and outdegree exponent γout = 2.87. To avoid
problems with controllability, extra edges are added randomly among different
strongly connected components, until strong connectivity is achieved on the entire
graph. The difference in the in/out exponent is still clearly visible. (b): Histogram
of outdegree/indegree ratio for the various networks. The SF networks reflect our
choice of γout < γin .

41

×10 10
0.1

8

tr(Wm )/n

λmin (Wm )

0.08
0.06
0.04
0.02

4

2

200

300

400

500

600

λmax(Wm )/λmin (Wm )

100

tr(W-1
)/n
m

6

10 10

10 5

10 0
100

200

300

400

500

100

200

100

200

300

400

500

600

300

400

500

600

10 25

10 20

10 15

600

n. inputs

n. inputs

10 1

ratio of tr(W m )

700
600
500
400
300
200
100
0

200

400

ratio of tr(W -1
)
m

10 0

10 -5

0

200

400

10 0

600

600

ratio of λmax(Wm )λmin (Wm )

ratio of λmin (Wm )

(a)

n. inputs

0

200

400

600

0

200

400

600

10 0

10 -1

10 -2

n. inputs

(b)

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 rw = wout /win (red) and
the same measure for random driver node assignments (blue). As the ratios in
(b) show, the improvement in all measures is normally of several orders of magnitude. Also the condition number of Wm (i.e., λmax (Wm )/λmin (Wm )) improves
substantially.

42

tr(Wz )

8

tr(W-1
)
z

6

random

10 8

w out /w in

4

tr(W-1
)/n
z

λ min (Wz )

10

tr(Wz)/n

λ min (Wz)

×10 -4
12

10 -1

10

6

10 4

2
100

200

300

400

500

600

100

n. inputs

200

300

400

500

600

100

n. inputs

200

300

400

500

600

n. inputs

(a)
1

overlap with w out /w in

0.6

overlap with w in/w out

1

overlap with λ min (Wz )

%

%

%

0.4
0.5

0.5

0.2
0

0
0

200

400

n. inputs

600

0
0

200

400

n. inputs

600

0

200

400

600

n. inputs

(b)

Figure S6: Driver node placement for a network of n = 1000 coupled harmonic
oscillators. The figure is the analogous of Fig. 4, but now the coupling matrix
K is fully connected. (a): Shown are means over 50 realizations (with 100 edge
weight samples taken for each realization). Red: driver node placement based
on λmin (Wz ). Violet: placement based on tr(Wz ). Green: placement based on
tr(Wz−1 ). Cyan: placement based on wout /win . 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 (Wz ) and
tr(Wz−1 ) tend to behave similarly and so do wout /win and tr(Wz ) (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 wout /win and
tr(Wz ) (> 90%), while λmin (Wz ) and tr(Wz−1 ) correspond to different node ranking
patterns. None of the strategies orders nodes according to win /wout , as expected.

43

10

3

10

10 2
10

random placement

4

10

tr(W)/n

10

w out /w in placement

tr(W)/n

tr(W)/n

λ min (W) placement
4

10 3
10 2
10

40

60

80

100

4

10 3
10 2
10

40

n. inputs

60

80

100

40

n. inputs

60

80

100

n. inputs

(a)
λ min (W) placement

random placement

10 10

40

60

n. inputs

80

100

10 25

tr(W-1 )/n

tr(W-1 )/n

10 20

-1

tr(W )/n

10

w out /w in placement

20

10 10

40

60

n. inputs

80

100

10 20

40

60

80

100

n. inputs

(b)

Figure S7: Minimum energy control of power grids with varying damping coefficients. North EU power grid. This Figure complements Fig. 5(b) of the paper.
(a): Control energy for the metric tr(Wr ) when the driver nodes are placed according to λmin (Wr ) (left panel), wout /win (mid panel), or randomly (right panel).
The color code is a function of the damping coefficients, with the same convention
as in Fig. 5(a) of the paper. The values of tr(Wr ) are shown in solid lines, while
in dotted lines the values of tr(Wz ) are shown (suitably normalized to eliminate
the explicit dependence from tf ). Values are averages over 100 realizations. (b):
Control energy for the metric tr(Wr−1 ), with the same conventions as in (a).

44

w

λ min (W) placement

out

/w

in

placement

×10 -5

random placement

0.07
0.06

10 -10

0.02
50

100

0.01

0.05

λ min (W)

0.03

λ min (W)

0.04

λ min (W)

2.5

0.04 10 -10
0.03
50

0.02

100

60

80

100

40

120

1

50

100

0.5

0.01
40

2 10 -10
1.5

60

80

100

40

120

60

80

100

120

n. inputs

n. inputs

n. inputs

(a)
w out /w in placement
10 4

3

10 3

10

10 2

random placement
10

tr(W)/n

10 4

tr(W)/n

tr(W)/n

λ min (W) placement

10 2

4

10 3
10 2

10

10

10
40

60

80

100

120

40

n. inputs

60

80

100

120

40

n. inputs

60

80

100

120

n. inputs

(b)
λ min (W) placement

w out /w in placement

10 10

random placement

20

tr(W-1 )/n

10

tr(W-1 )/n

-1

tr(W )/n

10 20

10 10

10 15

10 0

10 0
40

60

80

n. inputs

100

120

10 20

40

60

80

n. inputs

100

120

40

60

80

100

120

n. inputs

(c)

Figure S8: Minimum energy control of power grids with varying damping coefficients. IEEE 300 bus test power network. (a): Control energy for the metric
λmin (Wr ) when the driver nodes are placed according to λmin (Wr ) (left panel),
wout /win (mid panel), or randomly (right panel). The color code is a function of
the damping coefficients, using the same convention as in Fig. 5(a) of the paper.
The values of λmin (Wr ) are shown in solid lines, while in dotted lines the values
of λmin (Wz ) are shown (suitably normalized to eliminate the explicit dependence
from tf ). Values are averages over 100 realizations. Data are missing when the
Gramian Wr is too close to singular in too many trials. (b): Control energy for
the metric tr(Wr ) when the driver nodes are placed according to λmin (Wr ) (left
panel), wout /win (mid panel), or randomly (right panel). The values of tr(Wr ) are
shown in solid lines, while in dotted lines the values of tr(Wz ) are shown. (c):
Control energy for the metric tr(Wr−1 ), with the same conventions as in (a) and
(b).

45

3.5

w out /w in placement

λ min (W) placement

10 -3

12

(W)

10 -8
-13

400

1

8

min

2
1.5 10

10

600

800

10

-4

-10

10 -14

6

λ

λ min (W)

10

14

3
2.5

×10 -11random

×10 -5

λ min (W)

×10 -3

400

4

0.5

600

800

400

500

600

700

800

300

900

6

10

5

10 -15

4
3

10 -20
400

2

600

800

1

2

300

placement

-11

400

500

600

700

800

900

300

400

500

600

700

800

900

n. inputs

n. inputs

n. inputs

(a)
λ min (W) placement

random placement

10 4

10 3

10 2

10

tr(W)/n

4

tr(W)/n

tr(W)/n

10

w out /w in placement

10 3

400

500

600

700

800

900

10 3

10 2

10 2

300

4

300

400

n. inputs

500

600

700

800

900

300

400

n. inputs

500

600

700

800

900

n. inputs

(b)
λ min (W) placement
15

10 5

400

500

600

700

n. inputs

800

900

tr(W-1 )/n

10 10

300

random placement
10 25

10 20

tr(W-1 )/n

tr(W-1 )/n

10

w out /w in placement

10 10

300

400

500

600

700

n. inputs

800

900

10 20

300

400

500

600

700

800

900

n. inputs

(c)

Figure S9: Minimum energy control of power grids with varying damping coefficients. French high/mid voltage power grid. (a): Control energy for the metric
λmin (Wr ) when the driver nodes are placed according to λmin (Wr ) (left panel),
wout /win (mid panel), or randomly (right panel). The color code is a function of
the damping coefficients, using the same convention as in Fig. 5(a) of the paper.
The values of λmin (Wr ) are shown in solid lines, while in dotted lines the values
of λmin (Wz ) are shown (suitably normalized to eliminate the explicit dependence
from tf ). Data are missing when the Gramian Wr is too close to singular (mostly
when driver nodes are chosen randomly, right column). (b): Control energy for
the metric tr(Wr ) when the driver nodes are placed according to λmin (Wr ) (left
panel), wout /win (mid panel), or randomly (right panel). The values of tr(Wr ) are
shown in solid lines, while in dotted lines the values of tr(Wz ) are shown. (c):
Control energy for the metric tr(Wr−1 ), with the same conventions as in (a) and
(b).

46

λ

×10 -10 min

w

(W) placement

×10 -11 out

/w in placement

×10 -15random

3 10 -10

600

800

5

(W)
min

10 -20
400

600

2.5

800

0.5

10 -20

2
1.5

λ

(W)
min
400

10

λ

λ min (W)

2
1.5
1

3

15

2.5

placement

400

1

600

800

0.5

400

500

600

700

400

800

500

600

700

800

400

n. inputs

n. inputs

500

600

700

800

n. inputs

(a)

10 3

10 3

10 2

w out /w in placement

random placement
10 4

tr(W)/n

10 4

tr(W)/n

tr(W)/n

λ min (W) placement
10 4

10 2
10

10
400

500

600

700

800

10 3
10 2
10

400

500

n. inputs

600

700

800

400

500

n. inputs

600

700

800

n. inputs

(b)
λ min (W) placement

w out /w in placement

random placement

10

10 30

tr(W-1 )/n

10 20

tr(W-1 )/n

tr(W-1 )/n

10 30

10 25

10 25

10

400

500

600

n. inputs

700

800

400

500

600

n. inputs

700

800

400

500

600

700

800

n. inputs

(c)

Figure S10: Minimum energy control of power grids with varying damping coefficients. USA power grid. (a): Control energy for the metric λmin (Wr ) when the
driver nodes are placed according to λmin (Wr ) (left panel), wout /win (mid panel),
or randomly (right panel). The color code is a function of the damping coefficients,
using the same convention as in Fig. 5(a) of the paper. The values of λmin (Wr ) are
shown in solid lines, while in dotted lines the values of λmin (Wz ) are shown (suitably normalized to eliminate the explicit dependence from tf ). Data are missing
when the Gramian Wr is numerically too close to singular in too many trials. (b):
Control energy for the metric tr(Wr ) when the driver nodes are placed according
to λmin (Wr ) (left panel), wout /win (mid panel), or randomly (right panel). The
values of tr(Wr ) are shown in solid lines, while in dotted lines the values of tr(Wz )
are shown. (c): Control energy for the metric tr(Wr−1 ), with the same conventions
as in (a) and (b).
47

