Data-Driven Control of Complex Networks
Giacomo Baggio Danielle S. Bassett Fabio Pasqualetti
eess.SY cs.SY math.OC physics.soc-ph

naturemag


Department of Information Engineering, University of Padova, Padova, Italy

Departments of Bioengineering, Physics & Astronomy, Electrical
  & Systems Engineering, 
  Neurology, and Psychiatry, University of Pennsylvania, Philadelphia, USA
 Santa Fe Institute, Santa Fe, USA

Department of Mechanical Engineering, University of
  California at Riverside, Riverside, USA
 To whom correspondence should be addressed: mailto:fabiopas@engr.ucr.edufabiopas@engr.ucr.edu




 





Our ability to manipulate the behavior of complex networks depends
  on the design of efficient control algorithms and, critically, on
  the availability of an accurate and tractable model of the network
  dynamics. While the design of control algorithms for network systems
  has seen notable advances in the past few years, knowledge of the
  network dynamics is a ubiquitous assumption that is difficult to
  satisfy in practice, especially when the network topology is large
  and, possibly, time-varying. In this paper we overcome this
  limitation, and develop a data-driven framework to control a complex
  dynamical network optimally and without requiring any knowledge of
  the network dynamics. Our optimal controls are constructed using a
  finite set of experimental data, where the unknown complex network
  is stimulated with arbitrary and possibly random inputs. In addition
  to optimality, we show that our data-driven formulas enjoy favorable
  computational and numerical properties even compared to their
  model-based counterpart. Although our controls are provably correct
  for networks with linear dynamics, we also characterize their
  performance against noisy experimental data and in the presence of
  nonlinear dynamics, as they arise when mitigating cascading failures in power-grid networks and when manipulating neural activity in brain networks.



Data-Driven Control of Complex Networks
    Fabio Pasqualetti
    December 30, 2023
=======================================









§ INTRODUCTION

With the development of sensing, processing, and storing capabilities of modern sensors, massive volumes of information-rich data are now rapidly expanding in many physical and engineering domains, ranging from robotics <cit.>, to biological <cit.> and economic sciences <cit.>. Data are often dynamically generated by complex interconnected processes, and encode key information about the structure and operation of these networked phenomena. Examples include temporal recordings of functional activity in the human brain <cit.>, phasor measurements of currents and voltages in the power distribution grid <cit.>, and streams of traffic data in urban transportation networks <cit.>. When first-principle models are not conceivable, costly, or difficult to obtain, this unprecedented availability of data offers a great opportunity for scientists and practitioners to better understand, predict, and, ultimately, control the behavior of real-world complex networks.

Existing works on the controllability of complex networks have focused exclusively on a model-based setting <cit.>, although, in practice, constructing accurate models of large-scale networks is a challenging, often unfeasible, task <cit.>. In fact, errors in the network model (i.e., missing or extra links, incorrect link weights) are unavoidable, especially if the network is identified from data (see, e.g., <cit.> and Fig. <ref>(a)).  This uncertainty is particularly important for network controllability, since, as exemplified in Fig. <ref>(b)-(c), the computation of model-based network controls tends to be unreliable and highly sensitive to model uncertainties, even for moderate size networks <cit.>. It is therefore natural to ask whether network controls can be learned directly from data, and, if so, how well these data-driven control policies perform.


 


Data-driven control of dynamical systems has attracted increasing interest over the last few years, triggered by recent advances and
successes in machine learning and artificial intelligence <cit.>. The classic (indirect) approach to learn controls from data is to use a sequential system identification and control design procedure. That is, one first identifies a model of the system from the available data, and then computes the desired controls using the estimated model <cit.>. However, identification algorithms are sometimes inaccurate and time-consuming, and several direct data-driven methods have been proposed to bypass the identification step <cit.>. These include, among others, (model-free) reinforcement learning <cit.>, iterative learning control <cit.>, adaptive and self-tuning control <cit.>, and behavior-based methods <cit.>.

The above techniques differ in the data generation procedure, class of system dynamics considered, and control objectives. In classic reinforcement learning settings, data are generated online and updated under the guidance of a policy evaluator or reward
simulator, which in many applications is represented by an offline-trained (deep) neural network <cit.>. Iterative learning control is used to refine and optimize repetitive control tasks: data are recorded online during the execution of a task repeated multiple times, and employed to improve tracking accuracy from trial to trial. In adaptive control, the structure of the controller is fixed and a few control parameters are optimized using data collected on the fly. A widely known example is the auto-tuning of PID controllers <cit.>. Behavior-based techniques exploit a trajectory-based (or behavioral) representation of the system, and data that typically consist of a single, noiseless, and sufficiently long input-output system trajectory <cit.>. Each of the above data-driven approaches has its own limitations and merits, which strongly depend on the intended application area. However, a common feature of all these approaches is that they are tailored to or have been employed for closed-loop control tasks, such as stabilization or tracking, and not for finite-time point-to-point control tasks.


In this paper, we address the problem of learning from data point-to-point optimal controls for complex dynamical networks. Precisely, following recent literature on the controllability of complex networks <cit.>, we focus on control policies that optimally steer the state of (a subset of) network nodes from a given initial value to a desired final one within a finite time horizon. To derive analytic, interpretable results that capture the role of the network structure, we consider networks governed by linear dynamics, quadratic cost functions, and data consisting of a set of control experiments recorded offline. Importantly, experimental data are not required to be optimal, and can even be generated through random control experiments. In this setting, we establish closed-form expressions of optimal data-driven control policies to reach a desired target state and, in the case of noiseless data, characterize the minimum number of experiments needed to exactly reconstruct optimal control inputs. Further, we introduce suboptimal yet computationally simple data-driven expressions, and discuss the numerical and computational
advantages of using our data-driven approach when compared to the classic model-based one. Finally, we illustrate with different numerical studies how our framework can be applied to restore the correct operation of power-grid networks after a fault, and to characterize the controllability properties of functional brain networks.

While the focus of this paper is on designing optimal control inputs, the expressions derived in this work also provide an
alternative, computationally reliable, and efficient way of analyzing the controllability properties of large network systems. This
constitutes a significant contribution to the extensive literature on the model-based analysis of network controllability, where the
limitations imposed by commonly used Gramian-based techniques limit the investigation to small and well-structured network topologies <cit.>.




§ RESULTS




 §.§ Network dynamics and optimal point-to-point control
 We consider networks governed by linear time-invariant dynamics


    x(t+1)    = A x(t) + B u(t),  
    
    	y(t)    = C x(t),


where x(t)∈^n, u(t)∈^m, and y(t)∈^p denote, respectively, the state, input, and output of the network at
time t. The matrix A∈^n× n describes the (directed and weighted) adjacency matrix of the network, and the  matrices B∈^n× m and C∈^p× n, respectively, are typically chosen to single out prescribed sets of input and output nodes of the network. 

In this work, we are interested in designing open-loop control policies that steer the network output y(t) from an initial value
y(0)=y_0 to a desired one y(T)=yf in T steps. If yf is output controllable <cit.> (a standing assumption in this paper), then the latter problem admits a solution and, in fact, there are many ways to accomplish such a control task. Here, we assume that the network is initially relaxed (x(0) = 0), and we seek the control input
u^⋆_0:T-1 = [u^⋆(T-1)^⋯   u^⋆(0)^]^ that drives the output of
the network to yf in T steps and, at the same time, minimizes a prescribed quadratic combination of the control effort and locality of the controlled trajectories.

Mathematically, we study and solve the following constrained minimization problem: 
 

    u^⋆_0:T-1 = min_u_0:T-1       y_1:T-1^  Q  y_1:T-1 + u_0:T-1^  R  u_0:T-1
         s.t.  (<ref>)  and   y_T = yf,


where Q≽ 0 and R≻ 0 are tunable matrices[We let A≻ (≽)  0 denote a positive definite (semi-definite) matrix, and A^ the transpose of A.] that penalize output deviation and input usage, respectively, and subscript ·_t_1:t_2 denotes the vector containing the samples of a trajectory in the time window [t_1,t_2], t_1≤ t_2 (if t_1=t_2, we simply write ·_t_1). If Q=0 and R=I, then u^⋆_0:T-1 coincides with the minimum-energy control to reach yf in T steps <cit.>. 

Equation (<ref>) admits a closed-form solution whose computation requires the exact knowledge of the network matrix A and suffers from numerical instabilities (). In the following section, we address this limitation by deriving model-free and reliable expressions of u^⋆_0:T-1 that solely rely on experimental data collected during the network operation.








 §.§ Learning optimal controls from non-optimal data
 
We assume that the network matrix A is unknown and that N control experiments have been performed with the dynamical network in (<ref>). The i-th experiment consists of generating and applying the input sequence u_0:T-1^(i), and measuring the resulting output trajectory y_0:T^(i) (Fig. <ref>(a)). Here, as in, e.g., <cit.>, we consider episodic experiments where the network state is reset to zero before running a new trial, and refer to the  for an extension to the non-episodic setting and to the case of episodic experiments with non-zero initial state resets. We let U_0:T-1, Y_1:T-1, and Y_T denote the matrices containing, respectively, the experimental inputs, the output measurements in the time interval [1,T-1], and the output measurements at time T. Namely, 
 

    U_0:T-1   = [u_0:T-1^(1)  ⋯  u_0:T-1^(N)], 
    
    	Y_1:T-1   = [y_1:T-1^(1)  ⋯  y_1:T-1^(N)], 
     
    	Y_T   = [y^(1)_T  ⋯  y^(N)_T].


An important aspect of our analysis is that we do not require the input experiments to be optimal, in the sense of (<ref>), nor do we investigate the problem of experiment design, i.e., generating data that are “informative” for
our problem. In our setting, data are given, and these are generated from arbitrary, possibly random, or carefully chosen
experiments.

By relying on the data matrices in (<ref>), we derive the following data-driven solution to the minimization problem in (<ref>) (see the ):


    û_0:T-1 = U_0:T-1(I-K_Y_T(LK_Y_T)^†L)Y_T^† yf,


where L is any matrix satisfying L^L=Y_1:T-1^QY_1:T-1 + U_0:T-1^RU_0:T-1, K_Y_T denotes a matrix whose columns form a basis of the kernel of Y_T, and the superscript symbol ·^† stands for the Moore–Penrose pseudoinverse operation <cit.>.  
 


  §.§.§ Minimum number of data to learn optimal controls
 
Finite data suffices to exactly reconstruct the optimal control input via the data-driven expression in (<ref>) (see the ). In Fig. <ref>(c), we illustrate this fact for the class of Erdös–Rényi networks of Fig. <ref>(b). Specifically, the data-driven input û_0:T-1 equals the optimal one u^⋆_0:T-1 (for any target yf) if the data matrices in (<ref>) contain mT linearly independent experiments; that is, if U_0:T-1 is full row rank (Fig. <ref>(c), left). We stress that linear independence of the control experiments is a mild condition that is normally satisfied when the experiments are generated randomly. Further, if the number of independent trials is smaller than mT but larger than or equal to p, the data-driven control û_0:T-1 still correctly steers the network output to yf in T steps (Fig. <ref>(c), right), but with a cost that is typically larger than the optimal one. In this case, û_0:T-1 is a suboptimal solution to (<ref>), which becomes optimal (for any yf) if the collected data contain p independent trials that are optimal as well.
 
 




  §.§.§ Data-driven minimum-energy control
 By letting Q=0 and R=I in (<ref>), we recover a data-driven expression for the T-step minimum-energy control to reach yf. We remark that the family of minimum-energy controls has been extensively employed to characterize the fundamental capabilities and limitations of controlling networks, e.g., see <cit.>. After some algebraic manipulations, the data-driven minimum-energy control input can be compactly rewritten as (see the )


    û_0:T-1 
    	= (Y_TU_0:T-1^†)^† yf.


The latter expression relies on the final output measurements only (matrix Y_T) and, thus, it does not exploit the full output data
(matrix Y_1:T-1).  Equation (<ref>) can be further approximated as


    û_0:T-1≈ U_0:T-1 Y_T^† yf.


This is a simple, suboptimal data-based control sequence that correctly steers the network to yf in T steps, as long as p independent data are available. Further, and more importantly, when the input data samples are drawn randomly and independently from a distribution with zero mean and finite variance, (<ref>) converges to the minimum-energy control in the limit of infinite data (see the ). 

Fig. <ref>(a) compares the performance (in terms of control effort and error in the final state) of the two data-driven expressions (<ref>) and (<ref>), and the model-based control as a function of the data size N. While the data-driven control in (<ref>) becomes optimal for a finite number of data (precisely, for N=mT), the approximate expression (<ref>) tends to the optimal control only asymptotically in the number of data (Fig. <ref>(a), left). In both cases, the error in the final state goes to zero after collecting N=p data (Fig. <ref>(a), right). For the approximate control (<ref>), we also establish upper bounds on the size of the dataset to get a prescribed deviation from the optimal cost in the case of Gaussian noise. Our non-asymptotic analysis indicates that this deviation is proportional to the worst-case control energy required to reach a unit-norm target. This, in turn, implies that networks that are “easy” to control require fewer trials to attain a prescribed approximation error (see the ).




  §.§.§ Numerical and computational benefits of data-driven
  controls
 By relying on the same set of experimental data, in
Fig. <ref>(b), we compare the numerical accuracy, as
measured by the error in the final state, of the data-driven controls
 (<ref>) and (<ref>) and the
minimum-energy control computed via a standard two-step approach
comprising a network identification step followed by model-based
control design. First, we point out that if some nodes of the network are not accessible (C I) and no prior information about the network structure is available, then it is impossible to exactly reconstruct the network matrix A using (any number of) data <cit.>. In contrast, the computation of minimum-energy inputs is always feasible via our data-driven expression, provided that enough data are collected. We thus focus on the case in which all nodes can be
accessed (C=I). We consider Erdös–Rényi networks with n nodes as in
Fig. <ref>(b) and we select m= ⌊ n/10  ⌋ control nodes (forming matrix B). To reconstruct the network matrices A and
B, the subspace-based identification technique described in
. Data-driven strategies significantly outperform the standard
sequential approach for both dense (Fig. <ref>(b),
top) and sparse topologies (Fig. <ref>(b),
bottom). This poor performance of the standard approach is somehow expected because, independently of the
network identification procedure, the standard two-step approach
requires a number of operations larger than those required by the
data-driven approach, resulting in a potentially higher sensitivity to
round-off errors.  Also, it is interesting to note that the
data-driven approach is especially effective for large, dense networks
for which the standard approach leads to errors of considerable
magnitude (up to approximately 10^2).


A further advantage in using data-driven controls over model-based ones arises when dealing with massive networks featuring a small fraction of input and output nodes. Specifically, in Fig. <ref>(c) we plot the time needed to numerically compute the data-driven and model-based controls as a function of the size of the network. We focus on Erdös–Rényi networks as in Fig. <ref>(b) of dimension n≥ 1000 with ⌊ n/100 ⌋ input and output nodes and a control horizon T=50. The model-based control input requires the computation of the first T-1 powers of A (). The computation of the data-driven expressions (<ref>) and (<ref>) involves, instead, linear-algebraic operations on two matrices (U_0:T-1 and Y_T) that are typically smaller than A when n is very large  (precisely, when T<n/m and N<n). 

Thus, the computation of the control input via the data-driven
approach is normally faster than the classic model-based computation
(Fig. <ref>(c), left). In particular, the
data-driven control (<ref>), although
suboptimal, yields the most favorable performance due to its
particularly simple expression. Finally, we note that the error in the final state committed by the data-driven controls is always upper bounded by 10^-5 and thus it has a negligible effect on the control accuracy (Fig. <ref>(c), right).






  §.§.§ Data-driven controls with noisy data

The analysis so far has focused on noiseless data. A natural question is how the data-driven controls behave in the case of noisy data. If the noise is unknown but small in magnitude, then the established data-driven expressions will deviate slightly from the correct values (see the ). However, if some prior information on the noise is known, this information can be exploited to return more accurate control expressions. A particularly relevant case is when data are corrupted by additive i.i.d. noise with zero mean and known variance.[The different types of noise are assumed to be zero-mean to simplify the exposition. With slight modifications, non zero-mean noise could also be accommodated by our approach.]

Namely, the available data read as
 

    U_0:T-1   = U̅_0:T-1 +Δ_U, 
    
    	Y_1:T-1   = Y̅_1:T-1 + Δ_Y, 
     
    	Y_T   = Y̅_T + Δ_Y_T,


where U̅_0:T-1, Y̅_1:T-1, Y̅_T denote the ground truth values and Δ_U, Δ_Y, and Δ_Y_T are random matrices with i.i.d. entries with zero mean and variance σ_U^2, σ_Y^2, and σ_Y_T^2, respectively. Note that the noise terms Δ_Y and Δ_Y_T may also include the contribution of process noise acting on the network dynamics. In this setting, it can be shown that the data-driven control (<ref>) and the data-driven minimum-energy controls (<ref>) and (<ref>) are typically not consistent; that is, they do not converge to the true control inputs as the data size tends to infinity (see the  for a concrete example). However, by suitably modifying these expressions, it is possible to recover asymptotically correct data-driven formulas (). The key idea is to add correction terms that compensate for the noise variance arising from the pseudoinverse operations. In particular, the asymptotically correct version of the data-driven controls (<ref>) and (<ref>) read, respectively, as


    û'_0:T-1   =   
                              (Y_TU_0:T-1^(U_0:T-1U_0:T-1^-
                              Nσ_U^2I)^†)^†yf
                              , 
    û”_0:T-1   =    U_0:T-1 Y_T^(Y_TY_T^-Nσ^2_Y_TI)^† yf,


where we used the fact that X^†=X^(XX^)^† for any matrix X <cit.>, and Nσ_U^2I and Nσ_Y_T^2I represent the noise-dependent correction terms. Note, in particular, that if the noise corrupts the output data Y_T only, then (<ref>) coincides with the original data-driven control (<ref>), so that no correction is needed. Similarly, if the noise corrupts the input data U_T only, then (<ref>) coincides with the data-driven control (<ref>). 







 §.§ Applications




To demonstrate the potential relevance and applicability of the data-driven framework presented thus far, we present two applications of our data-driven control formulas. 



  §.§.§ Data-driven fault recovery in power-grid networks
 




We address the problem of restoring the normal operation of a
  power-grid network after the occurrence of a fault which
  desynchronizes part of the grid. If not mitigated in a timely manner, such
  desynchronization instabilities may trigger cascading failures that
  can ultimately cause major blackouts and grid disruptions
  <cit.>. In our case study,
  we consider a line fault in the New England power grid network comprising 39 nodes
  (29 load nodes and 10 generator nodes), as depicted in
  Fig. <ref>(a), and we compute an optimal point-to-point control from data to recover the correct operation of the grid. A similar problem is solved in <cit.> using a more sophisticated control strategy which requires knowledge of the network dynamics. As in <cit.>, we
  assume that the phase δ_i and the (angular) frequency
  ω_i of each generator i obey the swing equation dynamics
  with the parameters given in <cit.> (except for
  generator 1 whose phase and frequency are fixed to a constant, cf.
  ). Initially, each generator operates at a locally stable
  steady-state condition determined by the power flow equations. At
  time t=2s, a three-phase fault occurs in the transmission line
  connecting nodes 16 and 17. After 0.5s the fault is cleared;
  however the generators have lost synchrony and deviate from their
  steady-state values (Fig. <ref>(b)). To recover the normal
  behavior of the grid, 0.5s after the clearance of the fault, we
  apply a short, optimal control input to the frequency of the
  generators to steer the state (phase and frequency) of the
  generators back to its steady-state value. The input is computed
  from data via (<ref>) using N=4000 input/state experiments
  collected by locally perturbing the state of the generators around
  its normal operation point (see also ). We consider data
  sampled with period T_s=2.5× 10^-4s, and set the control
  horizon to T=400 time samples (corresponding to 0.1s), R=I,
  and Q=ε I with ε = 0.01 to enforce locality
  of the controlled trajectories. As shown in Fig. <ref>(c),
  the data-driven input drives the state of the generators to a point
  close enough to the starting synchronous solution (left,
  inset) so as to asymptotically recover the correct operation
  of the grid (right). Notably, as previously discussed,
  the computation of the control input requires only pre-collected
  data, is numerically efficient, and optimal (for the linearized dynamics).
  More generally, this numerical study
  suggests that the data-driven strategy (<ref>) could represent
  a simple, viable, and computationally-efficient approach to control
  complex non-linear networks around an operating point.









  §.§.§ Controlling functional brain networks via fMRI snapshots
 
We investigate the problem of generating prescribed patterns of activity in functional brain networks directly from task-based functional magnetic resonance imaging (task-fMRI) time series. Specifically, we examine a dataset of task-based fMRI experiments related to motor activity extracted from the Human Connectome Project (HCP) <cit.> (see Fig. <ref>(a)). In these experiments, participants are presented with visual cues that ask them to execute specific motor tasks; namely, tap their left or right fingers, squeeze their left or right toes, and move their tongue. We consider a set of m=6 input channels associated with different task-related stimuli; that is, the motor tasks' stimuli and the visual cue preceding them. As in <cit.>, we encode the input signals as binary time series taking the value of 1 when the corresponding task-related stimulus occurs and 0 otherwise. The output signals consist of minimally pre-processed blood-oxygen-level-dependent (BOLD) time series associated with the fMRI measurements at different regions of the brain (see also ). In our numerical study, we parcellated the brain into p=148 brain regions (74 regions per hemisphere) according to the Destrieux 2009 atlas <cit.>. Further, as a baseline for comparison, we approximate the dynamics of the functional network with a low-dimensional (n=20) linear model computed via the approach described in <cit.>, which has been shown to accurately capture the underlying network dynamics. 

In Fig. <ref>(b), we plot the inputs (top) and outputs (center) of one subject for the first sequence of five motor tasks. The bottom plot of the same figure shows the outputs obtained by approximating the network dynamics with the above-mentioned linear model. In Fig. <ref>(c), we compare the performance of the minimum-energy data-driven control in (<ref>) with the model-based one, assuming that the network obeys the dynamics of the approximate linear model. We choose a control horizon T=100, form the data matrices in (<ref>) by sliding a window of fixed size T over the available fMRI data, and consider a set of 20 orthogonal targets  {y_f,i}_i=1^20  corresponding to eigenvectors of the estimated T-step controllability Gramian (see  for further details). The top plot of Fig. <ref>(c) reports the error (normalized by the output dimension) in the final state of the two strategies, while the bottom plot shows the corresponding control energy (that is, the norm of the control input). In the plots, the targets are ordered from the most (y_f,1) to the least (y_f,20) controllable. The data-driven and the model-based inputs exhibit an almost identical behavior with reference to the most controllable targets. As we shift towards the least controllable targets, the data-driven strategy yields larger errors but, at the same time, requires less energy to be implemented, thus being potentially more feasible in practice. Importantly, since the underlying brain dynamics are not known, errors in the final state are computed using the identified linear dynamical model. It is thus expected that data-driven inputs yield larger errors in the final state than model-based inputs, although these errors may not correspond to control inaccuracies when applying the data-driven inputs to the actual brain dynamics. Ultimately, our numerical study suggests that the data-driven framework could represent a viable alternative to the classic model-based approach (e.g., see <cit.>) to infer controllability properties of brain networks, and (by suitably modulating the reconstructed inputs) enforce desired functional configurations in a non-invasive manner and without requiring real-time measurements.






§ DISCUSSION

In this paper we present a framework to control complex dynamical networks from data generated by non-optimal (and possibly random) experiments. We show that optimal point-to-point controls to reach a desired target state, including the widely used minimum-energy control input, can be determined exactly from data. We provide closed-form and approximate data-based expressions of these control inputs and characterize the minimum number of samples needed to compute them. Further, we show by means of numerical simulations that data-driven inputs are more accurate and computationally more efficient than model-based ones, and can be used to analyze and manipulate the controllability properties of real networks.

More generally, our framework and results suggest that many network control problems may be solved by simply relying on experimental data, thus promoting a new, exciting, and practical line of research in the field of complex networks. Because of the abundance of data in modern applications and the computationally appealing properties of data-driven controls, we expect that this new line of research will benefit a broad range of research communities, spanning from engineering to biology, which employ control-theoretic methods and tools to comprehend and manipulate complex networked phenomena.

Some limitations of this study should also be acknowledged and discussed. First, in our work we consider networks governed by linear dynamics. On the one hand, this is a restrictive assumption since many real-world networks are inherently nonlinear. On the other hand, linear models are used successfully to approximate the behavior of nonlinear dynamical networks around desired operating points, and capture more explicitly the impact of the network topology. Second, in many cases a closed-loop control strategy is preferable than a point-to-point one, especially if the control objective is to stabilize an equilibrium when external disturbances corrupt the dynamics. However, we stress that point-to-point controls, in addition to being able to steer the network to arbitrary configurations, are extensively used to characterize the fundamental control properties and limitations in networks of dynamical nodes. For instance, the expressions we provide for point-to-point control can also lead to novel methods to study the energetic limitations of controlling complex networks <cit.>, select sensors and actuators for optimized estimation and control <cit.>, and design optimized network structures <cit.>. Finally, although we provide data-driven expressions that compensate for the effect of noise in the limit of infinite data, we do not provide non-asymptotic guarantees on the reconstruction error. Overcoming these limitations represents a compelling direction of future work, which can strengthen the relevance and applicability of our data-driven control framework, and ultimately lead to viable control methods for complex networks.





§ METHODS




 §.§ Model-based expressions of optimal controls
 
The model-based solution to (<ref>) can be written in batch form as
 

    u^⋆_0:T-1 = (I - K_ C_T(MK_ C_T)^†M) C_T^†yf,


where C_T=[CB  CAB  ⋯ CA^T-1B] is the T-step output controllability matrix of the dynamical network in (<ref>), K_Y_T denotes a basis of the kernel of C_T, and M is any matrix satisfying M^M= H_T^Q H_T + R, with
 

    H_T = [       0       ⋯       ⋯       0      CB;       ⋮       ⋯       0      CB     CAB;       ⋮       ⋰       ⋰       ⋰       ⋮;       0      CB     CAB       ⋯ CA^T-2B;         ],


and 0 entries denoting p× m zero matrices. If Q=0 and R=I (minimum-energy control input), (<ref>) simplifies to u^⋆_0:T-1 =  C_T^†yf. Alternatively, if the network is target controllable, the minimum-energy input can be compactly written as
 

    u^⋆(t) = B^A^T-t-1C^ W_T^-1yf,    t=0,1,2,…,T-1.


where W_T denotes the T-step output controllability Gramian of the dynamical network in (<ref>) 
 

    W_T =  C_T C_T^ =∑_t=0^T-1 CA^tBB^(A^)^tC^,

which is invertible if and only if the network is target controllable.

Equation (<ref>) is the classic (Gramian-based) expression of the minimum-energy control input <cit.>. It is well-known that this expression is numerically unstable, even for moderate size systems, e.g., see <cit.>.




 §.§ Subspace-based system identification
 Given the data matrices U_0:T-1 and Y_T as defined in (<ref>) and assuming that C=I, a simple deterministic subspace-based procedure <cit.> to estimate the matrices A and B from the available data consists of the following two steps:


  * Compute an estimate of the T-step controllability matrix of the network as the solution of the minimization problem
 
  
    Ĉ_T = min_ C_TY_T- C_TU_0:T-1^2_F,

where ·_F denotes the Frobenius norm of a matrix. The solution to (<ref>) has the form
  Ĉ_T = Y_TU_0:T-1^†.


  * In view of the definition of the controllability matrix, obtain an estimate of the matrix B by extracting the first m columns of
  Ĉ_T. Namely,
  B̂ = [Ĉ_T,]_:,1:m, where [X]_:,i:j indicates the sub-matrix of X obtained from keeping the entries
  from the i-th to j-th columns and all of its rows. An estimate of the matrix A can be obtained as the solution to the least-squares problem
 

    Â = min_A [Ĉ_T]_:,m+1:mT - A  [Ĉ_T]_:,1:(T-1)m^2_F,


which yields the matrix Â=[Ĉ_T]_:,m+1:mT[Ĉ_T]_:,1:(T-1)m^†.

If the data are noiseless, the system is controllable in T-1 steps, and U_0:T-1 has full row rank, then this procedure provably returns correct estimates of A and B <cit.>.



 §.§ Power-grid network dynamics, parameters, and data generation
 

The short-term electromechanical behavior of generators {2,…,10} of the New England power-grid network are modeled by the swing equations <cit.>:

    δ̇_i   = ω_i,
    H_i/π f_bω̇_i    =  -D_iω_i + P_mi -G_iiE_i^2+∑_j=1,j i^10 E_iE_j(G_ijcos(δ_i-δ_j)+B_ijsin(δ_i-δ_j)).

where δ_i is the angular position or phase of the rotor in generator i  with respect to generator 1, and where ω_i  is the deviation of the rotor speed or frequency in generator i relative to the nominal angular frequency 2π f_b. The generator 1 is assumed to be connected to an infinite bus, and has constant phase and frequency. The parameters H_i and D_i are the  inertia constant and damping coefficient, respectively, of generator i. The parameter G_ii is the internal conductance of generator i, and G_ij + iB_ij (where i is the imaginary unit) is the transfer impedance between generators i and j. The parameter P_mi denotes the mechanical input power of generator i and E_i denotes the internal voltage of generator i. The values of parameters f_b, H_i, D_i, G_ij, B_ij, and P_mi in the non-faulty and faulty configuration are taken from  <cit.>, while the voltages E_i and initial conditions (δ_i (0), ω_i (0) = 0) are fixed using a power flow computation. In our numerical study, we discretize the dynamics (<ref>) using a forward Euler method with sampling time T_s=2.5× 10^-4s. Data are generated by applying a Gaussian i.i.d. perturbation with zero mean and variance 0.01 to each frequency ω_i. The initial condition of each experiment is computed by adding a Gaussian i.i.d. perturbation with zero mean and variance 0.01 to the steady-state values of δ_i and  ω_i of the swing dynamics (<ref>). 





 §.§ Task-fMRI dataset, pre-processing pipeline, and identification setup
 The motor task fMRI data used in our numerical study are extracted from the HCP S1200 release <cit.>. The details for data acquisition and experiment design can be found in <cit.>. The BOLD measurements have been pre-processed according to the minimal pipeline described in <cit.>, and, as in <cit.>, filtered with a band-pass filter to attenuate the frequencies outside the 0.06–0.12 Hz band. Further, as common practice, the effect of the physiological signals (cardiac, respiratory, and head motion signals) is removed from the BOLD measurements by means of the regression procedure in <cit.>. The data matrices in (<ref>) are generated via a sliding window of fixed length T=100 with initial time in the interval [-90,10]. We assume that the inputs and states are zero for times less than or equal to 10, i.e., the instant at which the first task condition is issued. We approximate the input-output dynamics with a linear model with state dimension n=20 computed using input-output data in the interval [0,150] and the identification procedure detailed in <cit.>.  When the estimated network matrix A has unstable eigenvalues, we stabilize A by diving it by ρ(A)+0.01, where ρ(A) denotes the spectral radius of A. Other identification parameters are as in <cit.>.




 §.§ Computational details
 
All numerical simulations have been performed via standard linear-algebra LAPACK routines available as built-in functions in
Matlab® R2019b, running on a 2.6 GHz Intel Core i5 processor with 8 GB of RAM. In particular, for the
computation of pseudoinverses we use the singular value decomposition method (command  in Matlab®) with a threshold of 10^-8.



 §.§ Materials and data availability
 
Data were provided (in part) by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University. The code and data used in this study are freely available in the public GitHub repository: https://github.com/baggiogi/data_driven_control.



§ ACKNOWLEDGEMENTS

This work was supported in part by awards AFOSR-FA9550-19-1-0235 and ARO-71603NSYIP, and by MIUR (Italian Minister for Education)  under the initiative “Departments of Excellence” (Law 232/2016).




§ AUTHOR CONTRIBUTIONS

G.B, D.S.B, and F.P. contributed to the conceptual and theoretical aspects of the study, wrote the manuscript and the . G.B. carried out the numerical studies and prepared the figures.



§ COMPETING INTERESTS

The authors declare no competing interests.






 


	
