Energy scaling of targeted optimal control of complex networks

Recently it has been shown that the control energy required to control a dynamical complex network is prohibitively large when there are only a few control inputs. Most methods to reduce the control energy have focused on where, in the network, to place additional control inputs. Here, in contrast, we show that by controlling the states of a subset of the nodes of a network, rather than the state of every node, while holding the number of control signals constant, the required energy to control a portion of the network can be reduced substantially. The energy requirements exponentially decay with the number of target nodes, suggesting that large networks can be controlled by a relatively small number of inputs as long as the target set is appropriately sized. We validate our conclusions in model and real networks to arrive at an energy scaling law to better design control objectives regardless of system size, energy restrictions, state restrictions, input node choices and target node choices.

From the Methods section we see that η may be computed from one target set size to another (which we call p min and p max ). To ensure that we compute a value of η that describes the entire network, we keep p min = 10% and compute values of logη p min →p max for larger values of p max . We see that the distributions as p max increases becomes 'sharper', i.e., that the standard deviation decreases, which is shown in the inset plot. After p max grows larger than 70%, we see that the improvement of the computed logη p min →p max slows down so that we do not need to compute η i for many additional points.  A three node network where node 1 is the driver, i.e., receives the control input, and node 3 is the target, i.e., the ouput of the system is the state of node 3, has an initial condition at the origin and a final condition when y f = x 3 (t f ) = 1. The solid lines correspond to minimum energy control, i.e., when Q 1 = Q 2 = O N andR = 1. The dashed lines correspond to a cost function where a weight of 1000 is included for the derivativesẋ 2 (t) andẋ 3 (t) and a small weight is included for the control input,R = 0.001. We can see that the rise of state one is more steep when a weight to the state derivative is included than for the minimum energy control trajectory. The state and control input weights can be tuned to achieve a desired state space trajectory. Legend -Targeting 100% * -Targeting 50% -Targeting 10% pβ β β . We show that for a variety of networks, real and model, scale-free and Erdos-Renyi, all nodes targeted or only some nodes targeted, the total energy for an arbitrary maneuverβ β β is well approximated by the open loop energy. Note that this approximation holds best when the controllability Gramian is poorly conditioned. Each control input is calculated for a cost function where Q 1 , Q 2 and R are appropriately dimensioned identity matrices. The model networks have the properties: n = 100, γ in = γ out = 3.0, k av = 5, and n d = 0.5. a Low average degree, k av = 2. b Moderate average degree k av = 5. The solid line has a slope of one. c Two real networks. The energy averaged over p control maneuvers β β β is shown in green for different values of p. The corresponding maximum energy is shown in red. Note that for any given p, the order of magnitude of the average energy is not much less than the order of magnitude of the maximum energy.
• -asc. in-degree -des. in-degree -asc. out-degree + -des. out-degree 1 Supplementary Figure 6. Effects of different selections strategies for the target nodes. We plot E (p) max versus the target fraction p/n for three real networks: the Florida foodweb [1], the Little Rock Lakes region foodweb [2], and a protein structure [3]. Nodes were removed from the target set in four different ways: (i) ascending in-degree, (ii) descending in-degree, (iii) ascending out-degree, (iv) descending out-degree. For each network n d = 0.45.  Both in the manuscript and here in the supplementary information, we examine how target control may benefit real networks compiled in datasets found throughout the scientific and engineering literature. We include the name, the reference, and some basic properties for each of the networks, as well as our computed value of η. In the table, n is the number of nodes, l is the number of edges, k av is the average degree, d is the diameter of the graph, and η is the scaling of the minimum control energy as we discuss in the manuscript and in the supplementary information.
Supplementary Note 1. Introduction to supplement 4 5 Complex networks have recently been used to model many distributed systems such as food webs, com-6 municating robots, financial interdependence, and social networks. While the dynamics of any one of those 7 networks are rich in nonlinearities and uncertain parameters, we will restrict ourselves to linear dynamics.

8
Linear dynamics are appropriate when a system is operating near a stable point, or if certain assumptions can 9 be made. Also, the differences between the specific dynamics make any overarching conclusions unlikely.

10
In the networks described before, often controlling every member is unnecessary which makes the control 11 action more 'expensive', by which we mean they require more effort, than is necessary. For instance, a preda-12 tor population in a foodweb may need to be reduced in order to improve a prey population, but other species 13 in the food web may not need be affected. In marketing, an ad agency may want to change the opinion of 14 a certain demographic, but not need to reach every member of the social network. A certain task, sent to a 15 robotic network may need to be performed by only a subset of its members. There are many control goals that 16 can be conceived of for complex networks where the desired final state should only be prescribed for some of 17 the members of the network but not for all of them, which we call target control. 18 We will show in the following sections that if target control is applicable to a dynamic network, the control 19 energy, or effort, decreases exponentially. We first provide a review of the minimum energy control problem 20 applied to a linear system with the addition of the concept of targeted states. Next, the exponential scaling  The fixed-end point minimum energy control problem is well-known in the optimal control field, especially 30 for a system described by linear dynamics, (1) What is less well known is the solution of the minimum energy control problem when the final condition is 32 only prescribed to some subset of the states. We introduce the minimum energy target control problem for 33 networks where the word target refers to those nodes with a prescribed final condition. The problem is as 34 follows: The matrix A ∈ R n×n is the adjacency matrix that describes the topology, or inter-connectedness, of the n 36 nodes, or states. The matrix B ∈ R n×m is the control input matrix that describes how the m control inputs are 37 distributed to the nodes. The matrix C ∈ R p×n is the output matrix that relates how each output is a linear 38 combination of the states. For the target control of complex networks formulation, we assume that B (C) has 39 columns (rows) that are all versors, i.e., each control input, u i (t), i = 1, . . . , m, is directed towards a single node 40 and each output, y j (t), j = 1, . . . , p, is the state of a single node (see Fig. 1a from the main manuscript for a 41 graphical description). The dynamical equation of an arbitrary node i is, where if there exists at least one coefficient b ik = 0 then node i is what we refer to as an input node. We will 43 assume that the system, (A, B,C), is output controllable so that, Each output is referred to as a targeted node. The solution of the minimization problem in Eq.
(2) is found 45 using Pontryagin's minimum principle [19] and is provided here both as a review and to establish how the 46 targeting aspect of our specific solution is applied. The Hamiltonian equation introduces n costates ν ν ν(t).
From the Hamiltonian equation, the following dynamical relations can be determined, The stationary equation is used to determine the optimal control input.
The time evolution of the costates can be determined in a straightforward manner with a final condition of the 50 form, ν ν ν(t f ) = C Tν ν ν f , whereν ν ν f ∈ R p as there are only p final conditions prescribed for the network.
With the optimal control input known, the time evolution of the states can also be determined, 52 The prescribed final condition for the targeted nodes is applied to determine the final, constant vectorν ν ν f .
If the system (A, B,C) is output controllable, then W is positive definite. When C has p rows (versors), the 55 matrix W p = CWC T , is the output controllability Gramian, and is a p × p principal submatrix of W .

56
Figures 2, 3, and 4 of the main text provide numerical evidence that the energy required for a control action 59 decreases exponentially as the number of target nodes decreases linearly. In the following derivation, we find 60 that the exponential decay of the energy is a result of a more fundamental property of the output controllability 61 Gramians W p . Here we show that for a broad class of networks and a random selection of the target nodes the 62 ratio of the smallest eigenvalues of two subsequent principal submatrices of the controllability Gramian W , by 63 which we mean the submatrices W p and W p−1 where W p−1 is W p after removing one additional row-column 64 pair, has a near constant value which we call η p = min{eig(W p−1 )}/ min{eig(W p )} ≈ constant. This is true for a typical sequence of random removals of target nodes (here by typical we mean that each node is assigned 66 the same probability of removal and the order of removal is random), while deviations from this behavior are 67 possible for specific removal strategies (see Section S6).

68
In the main text we have considered the average energy scaling when the cardinality of the target set de-69 creases from j to k, j > k. Here, we consider an iterative process as we remove one node at a time from the 70 target set. We say that two target node sets P p and P p+1 are adjacent if P p+1 = P p ∪ i and i / ∈ P p .

71
A symmetric, positive definite matrix W ∈ R n×n has principal submatrices W p ∈ R p×p , p < n where n − p 72 corresponding rows and columns of W have been removed. A principal submatrix, W p , has diagonal elements 73 which are also diagonal elements of the original matrix W . The eigenvalues of W p , µ (p) Consider the case when W p is W p+1 with one additional row-column pair removed, or in terms of the target sets, 76 P p ⊂ P p+1 which are adjacent. From Cauchy's interlacing theorem, the eigenvalues of W p thread between 77 the eigenvalues of W p+1 , The smallest eigenvalue of W p cannot be smaller than the smallest eigenvalue of W p+1 . We perform an iterative 79 process where at each step a row-column pair (without loss of generality here chosen to be the first row and 80 first column) is removed.
The matrixW p is a p × p principal submatrix of W p+1 with a first row of all zeros and a first column identical to 82 that of W p+1 . The matrix dW p consists of all zeros except for the first row which is identical to the first row of 83 W p+1 . The scalar w pp is the leading term in W p+1 and w p is the first column of W p+1 , after removing the entry 84 w pp . Note that the the set of eigenvalues ofW p is equal to the set of eigenvalues of W p with one additional 0 85 eigenvalue.

86
The smallest eigenvalue of W p+1 , µ (p+1) 1 , and the second smallest eigenvalue ofW p , µ 1 ( which is also the 87 smallest eigenvalue of W p ) are used to define the vectors v p+1 andv p , Pre-and post-multiplying Eq. (13) by v T p+1 andv p , respectively, will provide a relation between the smallest 89 eigenvalues of W p+1 and W p .
The matrix product W p+1 dW p is a matrix of all zeros except for the leading term which is one. Thus, the We use the definition of the 'worst-case' energy, E 1 to rewrite Eq. (16) in terms of energy, The last of Eq. (17) can be written in terms of any two target sets of size k and j, k < j and P k ⊂ P j , We defineη (k→ j) , which depends only on the two sets of target nodes P k and P j , as, In general, there are n! j!(n− j)! where i = p min , . . . , p max . It is seen through experiments that log η i is independent of the target set size (a 111 generic example is shown in Fig. 2) and can be computed for a given network. We stress that while we have i.e., not dependent on the current state. The closed-form solution is available after a similarity transformation 127 as explained below. The LQ optimal control problem is laid out below, The problem above considers systems with n states x i (t), i = 1, . . . , n, m control inputs u j (t), j = 1, . . . , m, 129 and p outputs y k (t), k = 1, . . . , p. We assign the same properties to B and C as we did previously, that the  H (x(t), ν ν ν(t), u(t)) = 1 2 The method defines three equations that, if satisfied, guarantees an optimal solution with respect to Eq. (25),

138
State Equation: The stationary equation provides the optimal control input, u u u * (t), What remains is to solve the dynamical system defined by the state and costate equations in Eq. (27). First,

140
Eq. (28) is applied to the state and costate equations to make the system homogeneous, The homogeneous Hamiltonian system is the following linear equation, state variables, while a final condition is imposed on the costate variables. Because of the coupled nature of 144 the problem, we cannot compute the state and costate trajectories individually. We use the relation between 145 the costate and the state, ν ν ν(t) = Sx(t) + ξ ξ ξ (t), where S is restricted to be symmetric, to define a similarity 146 transformation for the matrix in Eq. (31). We then obtain, The matrix in Eq. (32) has the following inversion property: We use the relation in Eq. (32) to rewrite Eq. (31) so thatξ ξ ξ (t) is decoupled from the states, x(t): The matrixÃ is defined as the augmented adjacency matrix because of its similar role in the state costate The matrixB acts as a control input matrix where ξ ξ ξ (t) acts as a pseudo control input, Our desire is to define S such thatQ is a zero matrix and so that the states, x(t), are decoupled from ξ ξ ξ (t), The solution for the alternate costate, ξ ξ ξ (t), is written in terms of a final condition, ξ ξ ξ (t f ) = ξ ξ ξ f .
All that remains is to apply the final output to the time evolution of the targets to define the vectorξ ξ ξ (t).
The matrixW = t f t 0 eÃ (t f −τ)B eÃ T (t f −τ) dτ is defined as the generalized controllability Gramian. We also can 161 define the control maneuver asβ β β = y f −CeÃ (t f −t 0 ) x 0 . The optimal control input is written in terms of the state 162 and alternate costate solutions, The optimal control input is the sum of a of the states and the final condition on the outputs, or more specifically, the targets.

166
The energy of the control action is defined as the cumulative effort of the control signal, u * c (t).
When the controllability Gramian is poorly conditioned, which is often the case in the control of large complex .
The vectorβ β β which we call the control maneuver, has information about the initial and final condition of the 173 system. We see from Eq. (43) that the optimal control input is actually the sum of two distinct components, u * c1 and u * c2 . However our approximation in Eq. (45) only takes into consideration u * c2 . We see that for both Comparing Eqs. (46) and (25) provides the relations Q = A T Q 1 A + Q 2 , M = A T Q 1 B, and R = B T Q 1 B +R. 182 We consider the formulation in Eq. (46) and provide an example in Fig. 3  at the end of its trajectory. More careful tuning can allow for many different shaped trajectories.

189
Here we will show that the order of the worst-case energy dominates the energy needed to reach any co-192 ordinate in state space. Consider a system with p targets, with initial condition at the origin and final output 193 located on the p-dimension unit hyper-sphere (this hyper-sphere is defined in the p-dimensional subspace of 194 phase space corresponding to the p target states). The average energy required to reach a location on the 195 p-dimensional unit hyper-sphere is determined from, where the vector v v v p ∈ R p and ||v v v p || = 1. Note that the vector β β β , which we call the control maneuver, can be 197 written as a weighted linear combination of the normalized eigenvectors of W −1 p , We use the summation in Eq. (48) to write Eq. (47) in terms of the eigenvectors of W −1 p .
For scale-free networks, we have seen the largest eigenvalue of W −1 p , 1 µ 1 , is far larger than any of the other So far we have investigated how the maximum energy is exponentially dependent on the target set size.

208
The relation log E (p) max = p n η + constant, where η is network specific, lets us compute the maximum energy if 209 η is known. To estimate η for a network, we have typically computed many values of η p , p = 1, . . . , n where 210 the p nodes in the target set are sampled randomly and the mean of η p is taken to beη. We now consider the 211 case if the target nodes are not chosen randomly, but instead, chosen with respect to the in-degree or out-degree 212 of the individual nodes. Figure 6 further investigates the dependece of E for the target nodes. Namely, the target nodes were chosen in order of ascending in-degree (AI), descend-214 ing in-degree (DI), ascending out-degre (AO), and descending out-degree (DO). As can be seen, when these 215 strategies are considered, E (p) max decreases in a way that strongly depends on the particular strategy applied and 216 substantially differs from network to network, i.e., it is network specific.