Optimal control of complex networks: Balancing accuracy and energy of the control action

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.¹⁵ 

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,²⁴ neuronal systems,^25,26 food webs,²⁷ and social systems.²⁸ 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.

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}, i=1,…,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 $κiin$ and the number of outgoing connections, called the out-degree $κiout$⁠. The average in-degree and average out-degree for a network are κ_av. 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., p(κ) ∼ κ^–γ where κ is either the in-degree or out-degree with corresponding γ_in and γ_out, and most often 2 ≤ γ ≤ 3.³¹ 

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.³² 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),…,xn(t)]T$ is the n × 1 time-varying state vector, $u(t)=[u1(t),…,um(t)]T$ is the m × 1 time-varying external control input vector, and $y(t)=[y1(t),…,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,³⁴ when all eigenvalues of A satisfy Reλ_i < 0 (See Section 3.3 in the book³⁴). 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⊆V$ 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).

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_n is the 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³⁵ (See sections 5.1 and 5.2 in the book³⁵). 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 with a final condition of the form, $ν(tf)=(1−α)CTν¯$⁠, where $ν¯=y(tf)−yf$

where $ν¯$ 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=∫t0tfeA(tf−τ)BBTeAT(tf−τ)dτ$ 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^T. The p dimensional vector $ν¯$ can be determined in a straightforward manner,

where $β=(CeA(tf−t0)x0−yf)$ and $Up=(α1−α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_p is positive semidefinite and the eigenvalues of U_p are the same as the eigenvalues of W_p plus the positive quantity $α1−α$⁠. Moreover, the eigenvectors of the matrices U_p and W_p are the same. From Eqs. (3)–(7), the optimal control input signal when the final condition is in the objective function is equal to

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.

The energy associated with the control input in Eq. (8), while only targeting the nodes for balanced control in $Pp$⁠, is defined as $ϵ(p)=∫t0tfu*(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−p)!$ 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 $β=CeA(tf−t0)x0−yf$ and $Mp=Up−1WpUp−1$ 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_p. Moreover, the following relation relates the eigenvalues of M_p and W_p:

where we denote the eigenvalues of M_p as $μi(p)$ and the eigenvalues of W_p as $λi(p), i=1,…,p$⁠. It follows that the energy expression (10) defines an ellipsoid in the variable $β$⁠. The axes of the ellipsoid are unaffected by the particular choice of the parameter α, while the width of each axis changes with the square root of the corresponding eigenvalue of the matrix M_p.

We consider the eigenvalues of M_p as $μi(p), i=1,…,p$⁠, which are ordered such that $0≤μ1(p)≤⋯≤μp(p)$⁠. By defining the magnitude of the vector, $|β|=β$⁠, we can define the ‘worst-case’ (or maximum) energy according to the Rayleigh quotient,

The upper extreme of the control energy denoted by $ϵmax(p)$⁠, for the control action in Eq. (8), is $max{ϵ(p)}∼μp(p)$⁠, which is what we call the ‘worst-case’ energy for optimal balanced control. For an arbitrary vector $β$⁠, 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, $μi(p)$⁠, which includes the worst-case energy.

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₀) = x₀ 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³⁰ 

The minimum energy for the optimal output cost control is

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_p, the output controllability Gramian, and the MBCE energy for this limit is

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₁ and target nodes 1, 2, and 3 are directly connected to outputs y₁, y₂, and y₃, 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⁻¹, the output states and the optimal control input are provided in panels (c) and (e), respectively. For small value of α = 10⁻⁷, 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 ϵ⁽³⁾≈ 4.2 and ϵ⁽³⁾≈ 219, respectively. We see that the energy can be reduced by 55 times in the former 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 $ζ=||y(tf)−yf||$⁠, the final state error at time t_f, $J1*=1−α2ζ2$⁠, and $J2*=α2ϵ(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 $−log α$ 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*→0, J2*/J*→1$ as α → 0. On the other hand, from Fig. 2, we see that as α approaches 1, $J2*/J*→0$ and $J1*/J*→1$⁠. As α increases (decreases), the error component (the energy component) becomes dominant in the optimal return value in Eq. (18).

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, 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_f = 1 and the fraction of input nodes n_d = 0.4. In Figs. 3 and 4, the solid line corresponds to the output cost control energy, 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 = γ_out). We notice that, when the network is more heterogeneous (γ 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 = γ_out = 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 (κ). 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_f for the same target nodes and in panels (e) and (f), we show that the optimal return function $J*$ decreases as α 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 = γ_out = 2.5, κ = 8. On the left half panels, we show the limiting relationship holds as the time horizon defined as t_f – t₀ changes. We notice that, when the time horizon is small (e.g., t_f – t₀ = 0.01), the terminal balancing is more beneficial compared to large time horizon (e.g., t_f – t₀ = 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_d changes. We notice that, when the number of designated input nodes is small (e.g., n_d = 0.4), the terminal balancing is more beneficial compared to a large number of input nodes (e.g., n_d = 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 ζ of the final state at final time t_f for the same target nodes and in panels (e) and (f), we show that the optimal return function $J*$ decreases as α 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,³⁷ Protein Structure,³⁷ Metabolic,³⁸ Food Web,^39–42 Social,^43–46 and Infrastructure.^47–49 For our simulation, we set t_f = 1 and n_d = 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 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 $ϵmax(p)$ versus α in Fig. 5. We see for small α, say, α = 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⁻¹, all of the networks need approximately the same amount of energy for balanced control.

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.³⁶ 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.

The controllability Gramian, W_p, can be calculated as a function of the eigendecomposition³⁰ of the state matrix A = VΛV⁻¹. For this article, we use the multi-precision package Advanpix for Matlab. The Matlab toolbox Advanpix⁵⁰ allows the computation of the eigendecomposition of W_p to be performed in an arbitrarily precise manner. This precision allows us to calculate the eigendecomposition of W_p, the invertible matrix U_p, and the matrix M_p accurately. We also use Advanpix when computing the energy for the cost function in Eq. (10).

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.

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

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.