Recently, it has been shown that the control energy required to control a large 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. We also have seen that by controlling the states of a subset of the nodes of a network, rather than the state of every node, 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. Here, we see that the control energy can be reduced even more if the prescribed final states are not satisfied strictly. We introduce a new control strategy called *balanced control* for which we set our objective function as a convex combination of two competitive terms: (i) the distance between the output final states at a given final time and given prescribed states and (ii) the total control energy expenditure over the given time period. We also see that the required energy for the optimal balanced control problem approximates the required energy for the optimal target control problem when the coefficient of the second term is very small. We validate our conclusions in model and real networks regardless of system size, energy restrictions, state restrictions, input node choices, and target node choices.

A typical problem in control applications is that of balancing the accuracy of the control action (i.e., how close one can reach a desired state) and the control effort (i.e., the expenditure of control signals). Here, we investigate the trade off between the control accuracy and the control effort in the context of dynamical complex networks.

## I. INTRODUCTION

Control of complex dynamical networks is not new.^{1–14} We define a network as controllable if a set of appropriate control signals can drive the network from an arbitrary initial condition to any final condition in finite time. If a network is controllable, there are several control signals we can choose from to approach the desired final condition. One important metric to characterize these control signals is the effort that each one requires. If **u**(*t*) is a control action which performs the task, the *control energy* associated with the control action is the cumulative magnitude of the control action over time. From optimal control theory, we can define the *optimal control action* which both satisfies our initial condition and final condition and minimizes the control effort required to perform the task over the time period of action.^{15}

Different types of control strategies in Refs. 6, 8–12, 14, and 16–19 have focused on controlling a broad range of networks such as power grids,^{20,21} communication networks,^{22,23} gene regulatory networks,^{24} neuronal systems,^{25,26} food webs,^{27} and social systems.^{28} While controlling a complex network, it becomes important to reduce the required control energy.

The minimum energy framework has been examined in Refs. 17 and 18 which have shown that based on the underlying network structure, the set of input nodes, the desired final state, and other parameters, the energy to control a network may lie on a distribution that spans a broad range of orders of magnitude.

One of the methods to reduce the required energy was investigated in Ref. 19, where additional control signals were added in optimal locations in the network according to each node's distance from the current set of control signals. Refs. 17 and 18 have only investigated the control energy for complex networks when the target set coincides with the set of all nodes. Refs. 14 and 29 examined methods to choose a minimal set of independent control signals necessary to control just the targets. Ref. 30, in which the control goal has been set to affect only a subset of the network nodes, chosen as the targets of the control action, has considered the effect of this choice on the required control energy. In fact, the energy requirements exponentially decay with the reduction of the number of target nodes.

With respect to the existing literature, here rather than exactly satisfying the prescribed final states we reduce the control energy by moving as close as possible to the prescribed final states. We introduce a new control strategy called *balanced control* where we set our performance measure as a convex combination of two objective functions: 1) a function which minimizes the distance from the final output states at final time to the desired final states and 2) a function which minimizes the control effort to achieve the first goal over the finite time period.

## II. RESULTS

### A. Problem formulation

Complex networks consist of two parts; a set of nodes with their interconnections that represent the topology of the network, and the dynamics which describe the time evolution of the network nodes. First, we summarize the definitions needed to describe a network. We define $V={i},\u2009i=1,\u2026,n$ to be the set of *n* nodes that constitute a network. The adjacency matrix is a real, square *n* × *n* matrix, *A*, which has nonzero elements *A _{ij}* if node

*i*receives a signal from node

*j*. For each node

*i*, we count the number of receiving connections, called the in-degree $\kappa iin$ and the number of outgoing connections, called the out-degree $\kappa iout$. The average in-degree and average out-degree for a network are

*κ*. One common way to characterize the topology of a network is by its degree distribution. Often the in-degree and out-degree distributions of networks that appear in science and engineering applications are scale-free, i.e.,

_{av}*p*(

*κ*) ∼

*κ*

^{–γ}where

*κ*is either the in-degree or out-degree with corresponding

*γ*and

_{in}*γ*, and most often 2 ≤

_{out}*γ*≤ 3.

^{31}

While most dynamical networks that arise in science and engineering are governed by nonlinear differential equations, the fundamental differences between individual networks and the uncertainty of precise dynamics make any substantial overarching conclusions difficult.^{6,29,31} Nonetheless, linear controllers have proven to be adequate in many applications by approximating nonlinear systems as linear systems in local regions of the *n*-dimensional state space.^{32} We examine linear dynamical systems, as it is a necessary first step to understand how balanced control may benefit nonlinear systems. The linear time invariant (LTI) network dynamics are

where $x(t)=[x1(t),\u2026,xn(t)]T$ is the *n* × 1 time-varying state vector, $u(t)=[u1(t),\u2026,um(t)]T$ is the *m* × 1 time-varying external control input vector, and $y(t)=[y1(t),\u2026,yp(t)]T$ is the *p* × 1 time-varying vector of outputs, or targets. The *n* × *n* matrix *A* = {*a _{ij}*} is the adjacency matrix described previously, the

*n*×

*m*matrix

*B*defines the nodes in which the

*m*control input signals are injected, and the

*p*×

*n*matrix

*C*expresses the relations between the states that are designated as the outputs. In addition, the diagonal values of

*A*,

*a*,

_{ii}*i*= 1,…,

*n*, which represent self-regulation, are chosen to be unique at each node (see justification in Ref. 33). These diagonal values are chosen to also guarantee that

*A*is Hurwitz so the system in Eq. (1) is internally stable. A matrix

*A*is called a Hurwitz,

^{34}when all eigenvalues of

*A*satisfy

**Re**

*λ*< 0 (See Section 3.3 in the book

_{i}^{34}). We restrict ourselves to the case when the matrix

*B*(

*C*) has linearly independent columns (rows) with a single nonzero element, i.e., each control signal is injected into a single node (defined as an input node) and each output is drawn from a single node (defined as a target node). We define $Pp\u2286V$ as the subset of target nodes and $p=|Pp|$ as the number of target nodes. A small sample schematic is shown in Fig. 1(a) that demonstrates the graphical layout of our problem emphasizing the graph structure and the role of input nodes (e.g., node 1 in the figure) and targets (e.g., node 1, 2, and 3 in the figure).

### B. Solution to the optimal balanced control problem

In our optimal balanced control problem, we attempt to minimize the following cost function,

subject to the constraints

Here, the final constraints are in the objective function and we call these constraints *soft constraints* as we do not require them to be satisfied exactly. Note that if we set *C* = *I _{n}*, where

*I*is the

_{n}*n*×

*n*identity matrix, then

**y**(

*t*) =

**x**(

*t*). The vector

**y**

_{f}is the prescribed final output state of the nodes described by the matrix

*C*. Here,

*α*∈ (0, 1] is a scaling parameter by which we can penalize the two performance measures in the cost function in (2) to balance the control energy. Note that in the case in which

*α*= 1, the cost function in Eq. (2) becomes the cost function associated with the optimal output cost control problem in Ref. 30, where different from Eq. (2), the final desired state is imposed as a

*hard constraint*.

The solution of the optimization problem in Eq. (2) is obtained using Pontryagin's maximum principle^{35} (See sections 5.1 and 5.2 in the book^{35}). The Hamiltonian equation introduces *n* costates $\nu (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 with a final condition of the form, $\nu (tf)=(1\u2212\alpha )CT\nu \xaf$, where $\nu \xaf=y(tf)\u2212yf$

where $\nu \xaf$ will be determined from the final output state. With the optimal control input known, the time evolution of the states can also be determined,

The final state of the targeted nodes can be determined

Here, $W=\u222bt0tfeA(tf\u2212\tau )BBTeAT(tf\u2212\tau )d\tau $ is the controllability Gramian. When *C* is defined as above, i.e., its rows are linearly independent versors, the *reduced Gramian W _{p}* is a

*p*-dimensional principal submatrix of

*W*, i.e., we write

*W*=

_{p}*CWC*. The

^{T}*p*dimensional vector $\nu \xaf$ can be determined in a straightforward manner,

where $\beta =(CeA(tf\u2212t0)x0\u2212yf)$ and $Up=(\alpha 1\u2212\alpha Ip+Wp)$. For 0 < *α* < 1, the *p* × *p* matrix *U _{p}* is always a symmetric, positive definite matrix, and invertible. In fact, the matrix

*W*is positive semidefinite and the eigenvalues of

_{p}*U*are the same as the eigenvalues of

_{p}*W*plus the positive quantity $\alpha 1\u2212\alpha $. Moreover, the eigenvectors of the matrices

_{p}*U*and

_{p}*W*are the same. From Eqs. (3)–(7), the optimal control input signal when the final condition is in the objective function is equal to

_{p}The equation for the time evolution of the outputs is equal to

Different from the formulation of linear optimal control commonly seen in texts on the subject, we approach the problem in two unique ways. First, we consider the control action that minimizes the cumulative magnitude of the control input, restricted by a left boundary condition applied to the state of the system (the initial condition) and a final condition applied to only the outputs of the system. Second, we make specific the methodology as it applies to networks by restricting our definitions of the matrices *B* and *C* as discussed previously.

### C. Optimal energy

The energy associated with the control input in Eq. (8), while only targeting the nodes for balanced control in $Pp$, is defined as $\u03f5(p)=\u222bt0tfu*(t)Tu*(t)dt$. Note that *ϵ*^{(p)} also depends on which *p* nodes are in the set, $Pp$. The energy *ϵ*^{(p)} is a measure of the ‘effort’ which must be provided to achieve the control goal. In the subsequent definitions and relations, when a variable is a function of *p*, we more specifically mean it is a function of a specific target set of size *p* of which there are $n!p!(n\u2212p)!$ possible sets. We can define the *minimum balanced control energy* (MBCE) when the control input is of the form in Eq. (8) as

where the vector $\beta =CeA(tf\u2212t0)x0\u2212yf$ and $Mp=Up\u22121WpUp\u22121$ is a *p* × *p* symmetric, real, semi-positive definite matrix. Note that the matrix *M _{p}* has the same set of eigenvectors as the matrix

*W*. Moreover, the following relation relates the eigenvalues of

_{p}*M*and

_{p}*W*:

_{p}where we denote the eigenvalues of *M _{p}* as $\mu i(p)$ and the eigenvalues of

*W*as $\lambda i(p),\u2009i=1,\u2026,p$. It follows that the energy expression (10) defines an ellipsoid in the variable $\beta $. The axes of the ellipsoid are unaffected by the particular choice of the parameter

_{p}*α*, while the width of each axis changes with the square root of the corresponding eigenvalue of the matrix

*M*.

_{p}### D. Worst case direction

We consider the eigenvalues of *M _{p}* as $\mu i(p),\u2009i=1,\u2026,p$, which are ordered such that $0\u2264\mu 1(p)\u2264\cdots \u2264\mu p(p)$. By defining the magnitude of the vector, $|\beta |=\beta $, we can define the ‘worst-case’ (or maximum) energy according to the Rayleigh quotient,

The upper extreme of the control energy denoted by $\u03f5max(p)$, for the control action in Eq. (8), is $max{\u03f5(p)}\u223c\mu p(p)$, which is what we call the ‘worst-case’ energy for optimal balanced control. For an arbitrary vector $\beta $, which can be represented as a linear combination of the eigenvectors of *M _{p}*, the energy can be defined as a weighted sum of the eigenvalues, $\mu i(p)$, which includes the worst-case energy.

### E. Solution to the optimal output cost control problem

Here, we summarize briefly the results for the optimal output cost control problem in Ref. 30. The cost function is minimized

where we choose **u**(*t*) such that it satisfies the prescribed initial state, **x**(*t*_{0}) = **x**_{0} and final output **y**(*t _{f}*) =

**y**

_{f}. We call this latter condition a hard constraint as these constraints need to be satisfied exactly.

The optimal control input signal when the final conditions are hard constraints (see the derivation in the supplementary material) and under the assumption that the triplet (A, B, C) is output controllable^{30}

The minimum energy for the optimal output cost control is

### F. Energy scaling with the penalizing factor *α*

We would like to determine the limiting energy when *α* → 0 in Eq. (10). At each *p*, $Pp$ contains *p* nodes in the target set and *U _{p}* is the

*p*×

*p*invertible matrix. When

*α*→ 0,

*U*→

_{p}*W*, the output controllability Gramian, and the MBCE energy for this limit is

_{p}We also provide the limiting behavior for the worst case energy direction

A small, three node example of the benefits of balanced control is shown in Fig. 1, where each node is in the target set (*p* = *n* = 3). Fig. 1(a) displays a sample network with three nodes. Input node (node 1) is in blue and target nodes for balanced control are in magenta (nodes 1, 2, and 3). Node 1 is directly connected to an input *u*_{1} and target nodes 1, 2, and 3 are directly connected to outputs *y*_{1}, *y*_{2}, and *y*_{3}, respectively. In panel (b), we examine the limiting relationship in Eq. (16) for the three node network. From Fig. 1(b), we also see how the balanced control strategy reduces the control energy as the penalizing factor *α* increases. For a large value of *α* = 10^{−1}, the output states and the optimal control input are provided in panels (c) and (e), respectively. For small value of *α* = 10^{−7}, the output states and the optimal control input are provided in panels (d) and (f), respectively. From panels (e) and (f), the integral of the magenta curves is *ϵ*^{(3) }≈ 4.2 and *ϵ*^{(3) }≈ 219, respectively. We see that the energy can be reduced by 55 times in the former case.

### G. Optimal return in limiting case

In the cost function (2), the two performance measures are multiplied each by a penalizing factor. It is important to investigate the relationship between the optimal return corresponding to each performance measure as *α* varies. The optimal return value corresponding to (2) can be written

where $\zeta =||y(tf)\u2212yf||$, the final state error at time *t _{f}*, $J1*=1\u2212\alpha 2\zeta 2$, and $J2*=\alpha 2\u03f5(p)$. We call the ratio $J1*/J*$ the optimal error return ratio and $J2*/J*$ the optimal energy return ratio. Note that the sum of the two ratios is equal to 1, i.e., $J1*/J*+J2*/J*=1$.

Figure 2 shows how the ratios $J1*/J*$ and $J2*/J*$ vary with *α* for the case of a scale free network with *n* = 300, *γ*_{in} = *γ*_{out} = 2.5, and *κ* = 8. We set the fraction of target nodes, *p*/*n* = 0.8 and the final time *t _{f}* = 1. The values on the abscissa axis are $\u2212log\u2009\alpha $ so that large values of

*α*are shown on the left hand side and small values of

*α*are shown on the right hand side. When

*α*→ 0,

*ζ*→ 0 faster than

*α*(which multiplies

*ϵ*

^{(p)}). Therefore, when

*α*is very small $J*$ is dominated by $J2*$, and we see from Fig. 2 that $J1*/J*\u21920,\u2009J2*/J*\u21921$ as

*α*→ 0. On the other hand, from Fig. 2, we see that as

*α*approaches 1, $J2*/J*\u21920$ and $J1*/J*\u21921$. As

*α*increases (decreases), the error component (the energy component) becomes dominant in the optimal return value in Eq. (18).

### H. Results and discussion

We perform numerical simulations to examine the two important results discussed in Subsections II E and II F. For our simulations in Figures 3 and 4, we consider scale-free model networks, constructed with the static model in Ref. 36 for specific parameters κ, the average degree, and *γ _{in}* =

*γ*=

_{out}*γ*, the power law exponent of the in- and out-degrees, and we choose the initial state at the origin,

**x**

_{0}=

**0**, and final state $||yf||=1$. We choose 10 different final states

**y**

_{f}uniformly distributed on the unit sphere centered in the origin, and we take the mean of all results over these realizations. We also consider

*t*= 1 and the fraction of input nodes

_{f}*n*= 0.4. In Figs. 3 and 4, the solid line corresponds to the output cost control energy,

_{d}*E*

^{(p)}, where

*p*/

*n*is the associated target set.

On the left half panels of Fig. 3, we construct scale free networks with the static model such that *n* = 300 and κ = 8 in each case. Each point is the average of 10 realizations and the bars represent one standard deviation. The target nodes are chosen from the set of nodes randomly and independently for each realization. In Fig. 3(a), the expected limiting relation discussed in Subsection II F is seen for each network irrespective of power-law exponent (*γ _{in}* =

*γ*). We notice that, when the network is more heterogeneous (

_{out}*γ*is low, e.g.,

*γ*= 2.5), the terminal balancing is more beneficial compared to the networks that are more homogeneous (

*γ*is high, e.g.,

*γ*= 3) as the balanced control energy remains close to the output control energy for homogeneous networks. We observe results qualitatively similar for any size of the terminal target set. On the left half panels of Fig. 3, we construct scale free networks with the static model such that

*n*= 300 and

*γ*=

_{in}*γ*= 2.5 in each case. Each point is the average of 10 realizations and the bars represent one standard deviation. Each set of target nodes is chosen from the set of nodes randomly and independently for each realization. In Fig. 3(b), the expected limiting relation discussed in Subsection II F is seen for each network irrespective of average degree (

_{out}*κ*). We notice that, when the network is more sparse (

*κ*is low, e.g.,

*κ*= 5), the terminal balancing is more beneficial compared to the networks that are dense (

*κ*is high, e.g.,

*κ*= 15) as the balanced control energy remains close to the output control energy for dense networks. This result holds for any size of the target sets. In panels (c) and (d), we show the error

*ζ*of the final state at final time

*t*for the same target nodes and in panels (e) and (f), we show that the optimal return function $J*$ decreases as

_{f}*α*decreases.

Besides the average degree and power-law exponent which describe the network (Fig. 3), there are other parameters that can affect the control energy such as the time horizon and the number of designated input nodes. In Fig. 4, each panel (a)–(f) corresponds to the size of target node set *p*/*n* = 0.8. The target nodes are chosen from the set of nodes randomly and independently for every realization. Each point is the average of 10 realizations and the error bars represent one standard deviation. The network has properties: *n* = 300, *γ _{in}* =

*γ*= 2.5, κ = 8. On the left half panels, we show the limiting relationship holds as the time horizon defined as

_{out}*t*–

_{f}*t*

_{0}changes. We notice that, when the time horizon is small (e.g.,

*t*–

_{f}*t*

_{0}= 0.01), the terminal balancing is more beneficial compared to large time horizon (e.g.,

*t*–

_{f}*t*

_{0}= 10) as the balanced control energy remains close to the output control energy. On the right half panels, we show the limiting relationship holds as the number of input nodes

*n*changes. We notice that, when the number of designated input nodes is small (e.g.,

_{d}*n*= 0.4), the terminal balancing is more beneficial compared to a large number of input nodes (e.g.,

_{d}*n*= 1) as the balanced control energy remains close to the output control energy for a large number of inputs. This result holds for any size of the target sets. In panels (c) and (d), we show the error

_{d}*ζ*of the final state at final time

*t*for the same target nodes and in panels (e) and (f), we show that the optimal return function $J*$ decreases as

_{f}*α*decreases.

We also analyze datasets collected from various fields in science and engineering to study how the worst-case energy for MBCE changes with *α* and the size of the target set under balanced control for networks with more realistic structures. In Fig. 5, we consider six groups of datasets: Circuit,^{37} Protein Structure,^{37} Metabolic,^{38} Food Web,^{39–42} Social,^{43–46} and Infrastructure.^{47–49} For our simulation, we set *t _{f}* = 1 and

*n*= 0.45. We take 30 realizations for one particular target fraction and take the mean over several realizations. We only show the results in Fig. 5 for

_{d}*p*/

*n*= 0.1 (for sufficient values of

*p*our results are qualitatively the same). For comparison among the real datasets, we choose one network from each network group and plot $\u03f5max(p)$ versus

*α*in Fig. 5. We see for small

*α*, say,

*α*= 10

^{−10}, the Metabolic network benefits most as the control energy for balanced control reduces significantly from the output control energy (magenta solid line). On the other hand, the Food Web and Social network are not benefited as much as balanced control energy remains approximately the same, in comparison. However, for large values of

*α*, say,

*α*= 10

^{−1}, all of the networks need approximately the same amount of energy for balanced control.

## III. METHODS

### A. Model networks

In our analyses, similar to Ref. 17, we assume that the networks have stable dynamics. The scale free model networks we consider throughout the paper are constructed with the static model.^{36} Diagonal noise, *δ _{i}*, is included, drawn from a uniform distribution between −1 and 1 so that the eigenvalues of the adjacency matrix are all unique. The weighted adjacency matrix

*A*is stabilized with a value

*η*such that each diagonal value of

*A*is {

*a*} =

_{ii}*δ*+

_{i}*η*where

*i*= 1,…,

*n*. The value

*η*is chosen such that the maximum eigenvalue of

*A*is equal to −1. The matrix

*B*is constructed by choosing which nodes in the network require an independent control signal. The matrices

*B*(

*C*) are composed of

*m*(

*p*) versors as columns (rows). We use the same method to choose input nodes as in Ref. 30.

### B. Numerical controllability

The controllability Gramian, *W _{p}*, can be calculated as a function of the eigendecomposition

^{30}of the state matrix

*A*=

*VΛV*

^{−1}. For this article, we use the multi-precision package Advanpix for Matlab. The Matlab toolbox Advanpix

^{50}allows the computation of the eigendecomposition of

*W*to be performed in an arbitrarily precise manner. This precision allows us to calculate the eigendecomposition of

_{p}*W*, the invertible matrix

_{p}*U*, and the matrix

_{p}*M*accurately. We also use Advanpix when computing the energy for the cost function in Eq. (10).

_{p}## IV. CONCLUSION

In this paper we provide an energy efficient control strategy we call balanced control strategy. We see that by changing the penalizing factor *α* in the cost function in Eq. (2), the control energy that is needed for balanced control can be reduced dramatically. For example, in Fig. 1(b), we see that the control energy can be reduced if we relax the final state conditions. We also see the limiting behavior as *α* → 0 MBCE approaches the output cost control energy. The above two results are general regardless of the network type, size, and other properties. See Figs. 1(b), 3, 4, and 5. From Figure 3, we see that the sparse and heterogeneous networks benefit more from our balanced control strategy than dense and homogeneous networks. We discuss the effect of other parameters, especially time horizon and number of input nodes on the MBCE and its limiting behavior. Several real datasets have also been examined to verify this result. In Fig. 5, we compare the results for different groups of real networks and conclude that the biological networks (e.g., metabolic, protein structure) are those that benefit most from the balanced control strategy.

## V. SUPPLEMENTARY MATERIAL

See supplementary material for complete derivation of the minimum energy output control problem and the definition of versor.

## ACKNOWLEDGMENTS

We gratefully acknowledge support from the National Science Foundation through NSF Grant No. CMMI- 1400193, NSF Grant No. CRISP- 1541148, and the Office of Naval Research through ONR Award No. N00014-16-1-2637.