A topology-dynamics-based control strategy for multi-dimensional complex networked dynamical systems

Complex systems are omnipresent and play a vital role in in our every-day lives. Adverse behavior of such systems has generated considerable interest in being able to control complex systems modeled as networks. Here, we propose a topology-dynamics-based approach for controlling complex systems modeled as networks of coupled multi-dimensional dynamical entities. For given dynamics and topology, we introduce an efficient scheme to identify in polynomial time a finite set of driver nodes, which – when endowed with the control function – steer the network to the desired behavior. We demonstrate the high suitability of our approach by controlling various networked multi-dimensional dynamics, coupled onto different topologies.

Most of these strategies solely depend on the system's network structure and have no direct connection to the underlying dynamics, and it only suffices for the dynamics to be describable by homogeneous dynamical equations.
The complexity of state-of-the-art strategies 32 for controlling linear time-invariant dynamics of networks with N nodes, and if one uses the Kalman rank condition, is 2 N times the complexity of evaluating the rank of the Kalman controllability matrix, which stems from the fact that there are 2 N − 1 possible combinations for selecting driver nodes. For strategies based on structural controllability, the minimum number of driver nodes N D , needed to fully control a directed network, is determined by the maximum matching in the network, where the unmatched nodes are exactly the driver nodes 8 . This allows one to find the driver nodes with the Hopcroft-Karp algorithm 33 , which has complexity O N L ( ) where L denotes the number of links. Strategies based on exact controllability can be applied to arbitrary network structures and link weights without any limitations, and with the Popov-Belevitch-Hautus rank condition 34 its complexity is O N N ( ln ( )) 2 2 . Another strategy, which also solely takes into account the structure of the underlying network, is based on the minimum dominating set 35,36 and has complexity . O(1 5137 ) N 37 . Here, we propose a control approach and a driver-node-identification scheme -based on a transformation of state variables of multi-dimensional dynamics -that makes use of both the equations of motion of the subsystems and topological properties of the underlying network. This enables us to determine a suitable control function. When endowed with this control function, the driver nodes steer the network to the desired state. Our approach belongs to the category (i) but our driver-node-identification scheme can also be combined with control strategies from category (iv) (Supplementary Information).

Results
We illustrate our approach -based on pinning -to control a complex dynamical network towards a desired state. Our approach relies on the linearization (as well as transformation of state variables) of any given dynamics in the vicinity of a desired final state and on meeting the necessary criteria to linearly stable state. In addition, we present an algorithm that allows one to determine -for the given dynamics and topology -a small set of the required driver nodes. Our approach is established upon few assumptions that will be stated in what follows and they are set to be default in the whole paper, unless explicitly stated. First, we assume that the desired state is a fixed point of the system. We note though that one can employ our approach also in cases where the desired state is a nominal trajectory by combining it with the master stability function (MSF) formalism 26,27 (Supplementary Information). Second, we assume full observability of the system, i.e., the state of the whole system is exactly known at any time t. Third, we assume accessibility of the system, i.e., each and any subsystem can be pinned by the specific control function.
control strategy. We present the basic ideas to determine the control function in the Methods section. Here, we illustrate our approach by controlling a network of second-order Kuramoto oscillators. The N phase oscillators θ i (i ∈ {1, …, N}), placed on a network, evolve according to 38 where ω i denotes the unique natural frequency (deviation from a nominal frequency) of the ith oscillator and α i is the damping constant (we fix α i = 2 for all oscillators and consider ω i to be uncorrelated random numbers, Gaussian distributed, with zero mean and fixed standard deviation). Also the non-zero coupling strengths λ ij are independent, normally distributed, random variables with a mean in the interval [0.1; 1.5] and a standard deviation amounting to one fifth of the mean. We have chosen the values of parameters such that they allow the system to attain both synchronised and asynchronous states with different initial conditions. By integrating Eq. (1) for a given network topology and initial condition, we derive the temporal evolution of variables θ i (t) and ν θ =  t t ( ) ( ) i i , which fully determine the system's dynamics. If the phase θ i (t) of oscillator i approaches the phase of the mean field (i.e., ψ = ∑ ∑ θ θ apart from a constant shift), it is locked to the mean field, and if this happens to all oscillators, the system is synchronised. In the left part of Fig. 1, we show an exemplary frequency dynamics of a network of non-identical second-order Kuramoto oscillators that, starting from some random initial condition, does not achieve the synchronised state. We note that for other initial conditions and network configurations the system may attain the synchronised state.
In the following, we refer to the synchronised and the incoherent state as the desired and the undesired state, respectively. Our approach consists of identifying the oscillators (or nodes) that cause positive eigenvalues of the network's Jacobian matrix of transformed state variables (see below, i.e. Eqs. 11 and 12), and controlling them will lead to a linearly stable desired state (Methods). To do so, we inspect the spectrum of eigenvalues around the desired state. If all eigenvalues have a negative real part, the system is linearly stable around that state. We verify the vanishing of the positive real part of eigenvalues by means of the Gershgorin circle theorem 15,39 (Methods). In order to guarantee that all components of the dynamical system will approach the desired state (i.e., all eigenvalues will lie in the negative half-plane of the complex plane), we move all Gershgorin disks to the left half-plane of the complex plane by adding a pinning term (in the general form of a feedback-control term) such as β where β i is the pinning strength and x i the state variable of node i. Here ⁎ x i denotes the i-th component of the fixed point. We note that, depending on the definition of state variables, the pinning term might find some nontrivial form as shown in two examples in the reminder of this paper.
The Gershgorin disks related to ν ∂ i of the corresponding Jacobian (Supplementary Information) have the following centers C i and corresponding radii R i : Accordingly, these problematic nodes or "driver nodes" need to be pinned.
Identification of driver nodes. The Gershgorin circle theorem provides us with a sufficient but not a necessary condition for controlling a network of oscillators. Applying the obtained control function to the set of Eq. (2) and using condition (4) can thus result in an abundant number of driver nodes. We therefore aim at identifying a finite and smaller number of driver nodes (set of size N D ≪ N) that mainly affect the network's dynamics. To do so, we propose the following scheme: • Step 1: add the control term with strength β 1 satisfying condition (4) to node/oscillator 1, calculate the maximum real part of the Jacobian's eigenvalues, and save its value as e 1 .  starting from some random initial condition from which the system does not achieve the fully synchronised state. (b) frequency dynamics of the controlled system with N D = 10 driver nodes starting from the same initial condition. The controlled system attains the synchronised state after about 20000 integration time steps (4thorder Runge-Kutta method with step size dt = 10 −3 ).

•
Step 3: identify the minimum value in this set. The corresponding node will be a driver node and its control term will be fixed. • Step 4: if min(e i ) ≥ 0, repeat steps 1-3 and identify other driver nodes.
Note that the linearization of the dynamics around a fixed-point may lead to a small error for finite deviations from the fixed-point. This error can be avoided by introducing a small negative upper bound for min(e i ). For the analysis presented below, we use the condition min(e i ) ≥ −0.2.
The computational complexity of our driver-node-identification scheme scales with network size as N δ × m 2.4 , where δ < 4.4 and m is the dimensionality of the state variable of a given node (e.g., m = 2 for the second-order Kuramoto oscillators). Applying our control strategy to the system presented in the upper part of Fig. 1, we find that only N D = 10 driver nodes (out of total 120 nodes) are required to steer the system dynamics from the undesired state (asynchronous) to the desired one (synchronised, cf. right part of Fig. 1).
impact of coupling topology. To demonstrate extendability of our observations beyond exemplary dynamics, we now investigate how the important aspects of our control strategy -namely the required number of driver nodes N D and its fraction n D = N D /N (pinning density) as well as the pinning strengths β -depend on properties of the coupling topology. We here consider Erdős-Rényi (ER) networks 40 (wiring probability p ∈ {0.03, …, 0.8}) and scale-free (SF) networks (Barabási-Albert model 41 with an average number of neighbors 4), both for a range of the number of nodes (N∈{50, …, 290}). We generate 10 realisations for each network of coupled second-order Kuramoto oscillators.
For ER networks, pinning density n D increases with network size N ( Fig. 2a), which can be related to the fact that the number of links in such networks increases with N as N 2 for a constant wiring probability p. Given that the number of driver nodes in our control strategy depends on the number of links, we expect N D /N 2 = n D /N vs. N to be constant (see the lower inset of Fig. 2a). We also find that n D saturates at about 30% for large wiring probabilities p (see the upper inset of Fig. 2a). Note that other final-state-oriented control strategies, for instance the feedback vertex set (FVS) control scheme 24 , require a two-to three-fold higher pinning density (from which one finds 75% for p = 0.03, 89% for p = 0.175, and 93% for p = 0.3).
For SF networks, pinning density n D is independent of the network size N and takes on values between 3% -6% (Fig. 2b). For such networks, the FVS control scheme requires 48% of nodes to be controlled, while it amounts to 37% with structural controllability (in Supplementary Information, we report also our findings obtained with the MSF formalism). The probability density of pinning strength β depends on the properties of the coupling topology (Fig. 2c). For SF networks, we find β to be only marginally influenced by the number of nodes N, and control of these networks requires comparably low strengths (β  2). For ER networks with for instance N = 120, we find an only marginal influence of the wiring probability (here p = 0.03, p = 0.3, and p = 0.8), but control of these networks requires higher pinning strengths β  20) in comparison with SF networks. characteristics of driver nodes. Having characterized important properties of our control approach (pinning density and pinning strength), we next investigate whether driver nodes distinguish themselves by specific characteristics. Our findings obtained from extensive numerical simulations and analyses of oscillatory dynamics (second-order Kuramoto oscillators) on about 1000 ER and SF networks are compiled into Fig. 3. Mean natural frequency of oscillators associated with driver nodes in ER networks equals the one associated with non-driver nodes (with ω 〈 〉  0) although with a reduced spread (~2.2). For SF networks, however, natural frequencies of driver node oscillators have a localized distribution around ~4 and ~−4, but the mean degrees of the driver nodes (〈k D 〉≃1.5) are smaller in comparison with the mean degree of all nodes (〈 〉  k 4). These findings demonstrate that properties of driver nodes -namely their natural frequency and their mean degree -are related to the underlying coupling topology. www.nature.com/scientificreports www.nature.com/scientificreports/ We also find that the pinning strength β is highly correlated to the degree k of driver nodes for both ER (with p = 0.03) and SF networks of second-order Kuramoto oscillators. The averaged values of degree k of driver nodes and pinning strength β for both networks are 〈〈 〉  k 11 and β 〈 〉  20.
controlling second-order kuramoto oscillators on a power grid. We now demonstrate the efficiency of our approach to control oscillators networks with realistic coupling topologies. To this end, we place second-order Kuramoto oscillators onto the relatively simple, but representative network of the IEEE Reliability Test System 42,43 (Fig. 4a, Supplementary Information), which is often used to investigate stability of real-world power grids, and contrast our findings to those obtained with the FVS control scheme. The second-order Kuramoto oscillators on this power grid will synchronise, if we choose a sufficiently high damping constant α (Eq. 1; Table S1 in Supplementary Information). Upon reducing the damping constant (from α = 2 to α = 0.8) and keeping other network properties unchanged, the grid will not synchronise (Fig. 4b). Using our driver node identification scheme, we find that only N D = 3 out of 73 nodes (4%) are required to control the grid, and when adding the control function with the form β ν i to the dynamics of these driver nodes, the grid converges to the fully synchronised state (Fig. 4c).
The FVS control scheme requires 53% of the grid's nodes to be controlled. In general, sets of driver nodes identified with either scheme do not fully overlap, and for the example presented above, only two driver nodes identified with our scheme are also driver nodes identified with the FVS scheme.
We note that the pinning strength β is related to parameters and characteristics of the power grid (Table S1, Supplementary Information) and for a given damping constant (dissipation coefficient) can be adjusted by modifying either the adjacency matrix or the reactance of connections linked to a driver node.
With an eye on applications to real power grids, we now aim at controlling the IEEE Reliability Test System under the impact of noise, an essential ingredient for power grids in view of increasing renewable energy sources [44][45][46][47][48][49][50][51][52] . Here, we consider the dynamics of the rotors (Supplementary Information, Eq. S.16) and add a noise term only to the generators (i.e., rotors with positive natural frequency) as follows www.nature.com/scientificreports www.nature.com/scientificreports/ where η i (t) is an uncorrelated white noise chosen randomly between ω i (1 − γ/2) and ω i (1 + γ/2) with uniform distribution ( Supplementary Information). Here, γ∈{0, 2} determines the standard deviation of the noise, which is temporally and spatially uncorrelated (w.r.t. the network).
The added noise η i (t) does not depend directly on the state variables of the system, and accordingly, does not appear in the Jacobian matrix, similar to the natural frequency ω i . Therefore, it can be considered as the momentary natural frequency ω i of generator i when added to the constant power input ω i . However, changing ω i at every time step changes the fixed point of the system as well as its corresponding Jacobian matrix. Since our driver-node-identification scheme is based on the eigenvalues of the Jacobian matrix at the given fixed point, one should run the proposed scheme at every time point. In practice, however, one can find the sets of driver nodes for a given number of noise realisations where at each realisation, the natural frequency of each node is drawn from the aforementioned distribution. Pinning the union of these sets of the driver nodes would be sufficient for stabilising the system when the pinning strength of every node is equal to the maximum of the pinning strengths for that node among all realisations. It should be noted that the noise in the natural frequency is directly related to the noise in the power input of the generators (see Supplementary Information for more details).
Here, we identify driver nodes for 500 realisations of the noise. For each realisation, the identified driver nodes coincide with the three nodes identified for the noise-free system. Nevertheless, the required pinning strengths (normalized to the respective values from the noise-free system) attain a wider distribution when the standard deviation of the noise is increased (Fig. 5), as expected. If we pin the union of sets of driver nodes for each realisation with the maximum pinning strength obtained from the different realisations, we can control the noisy system as shown in Fig. 5. controlling networks of complex multi-dimensional dynamical systems. Extending our results for relatively simple oscillators, we finally demonstrate the efficiency of our approach to control more complicated dynamical systems. As an example, we consider networks of N coupled Rössler oscillators, where oscillator i evolves according to x i denotes the oscillator's internal state and σ is the global coupling strength. We fix the model's control parameter (a = −0.2, b = 0.2, c = 5.4) so that the system either approaches its fixed point (here fully synchronised state) or diverges depending on initial conditions. Accordingly, the aim of our control scheme will be to steer the system to the fully synchronised state and to prevent divergence (in Supplementary Information, we present an alternative iterative scheme for identifying driver nodes if the desired state is a chaotic attractor).
Oscillators are coupled through their x 1 and x 3 components via the coupling function H, and g ij is the Laplacian matrix associated with the weighted adjacency matrix A of the network as , with a ij denoting elements of A.
To control this network, we derive a term (Supplementary Information) of the form (6)) -enables us to move the eigenvalues of the corresponding Jacobian to the left half-plane of the complex plane. Here ⁎ x j i (j∈{1, 2, 3}) are the stationary solutions of Eq. (6). Also, α in is the sum of weights connected to node i. The pinning strength β i of node i satisfies the following relation where the term J is related to the Jacobian of the controlled system (Supplementary Information). Note that the suggested control function depends on both, parameters controlling the dynamics (a and σ) and on properties of the underlying connection topology (k i in ). We now proceed as above to identify sets of driver nodes for ER and SF networks with N nodes and show in Fig. 6 how pinning density depends on network size N. In the case of ER networks, and in contrast to the networks of second-order Kuramoto oscillators, n D does not increase monotonically with increasing N, however, a maximum can be observed at N ≈ 130. Moreover, n D decreases monotonically with increasing the wiring probability p, which again contrasts our findings for the networks of second-order Kuramoto oscillators.

Figure 5.
Control of noisy second-order Kuramoto oscillators coupled onto a power grid. (a) Normalised pinning strengths β for the three identified driver nodes of the IEEE73 test system (see Fig. 4) under the impact of uncorrelated uniformly distributed white noise with different standard deviation as determined by γ (β 0 denotes the pinning strength of the noise-free system). Data for nodes 6 and 11 are slightly shifted to the right for the sake of visibility. (b) Noisy frequency dynamics (γ = 0.1) of second-order Kuramoto oscillators coupled onto the IEEE73 test system starting from some random initial condition from which it does not achieve the fully synchronised state. (c) Frequency dynamics of the controlled noisy power grid with N D = 3 driver nodes starting from the same initial condition. Pinning strengths of driver nodes amounted to β i ∈{53, 49, 103}, i∈{2, 6, 11}. A convergence of the power grid frequency is attained after about 50000 integration time steps (2ndorder stochastic Runge-Kutta method with step size dt = 10 −3 ). www.nature.com/scientificreports www.nature.com/scientificreports/ The existence of the maximum in the dependence of n D on N is related to the fact that the pinning density is affected by two variables -the mean degree and the network size. One can speculate that for small values of N, the impact of the network size dominates and leads to an increase of n D . Likewise, for larger values of N, the impact of the mean degree leads to a decrease of n D . To verify this hypothesis, one can study the dependence of n D vs N with a fixed mean degree. The Erdős-Rényi networks could not be used for this purpose since with a fixed wiring probability p, the mean degree increases with network size N. Accordingly, we use Watts-Strogatz (WS) networks (here with a rewiring probability of 0.9 and with 4 "nearest neighbors") which are random, sparse networks with a fixed mean degree. For such networks, n D indeed increases with N (lower inset in Fig. 6a). Therefore, by taking into account the decreasing behavior of n D with an increasing wiring probability p (and therefore increasing mean degree; upper inset in the Fig. 6a), and increasing of n D with N, we can understand the existence of the aforementioned maximum for ER networks. Finally, in the case of SF networks, we observe n D to increase monotonically with network size, in contrast to the SF networks of second-order Kuramoto oscillators.
Investigating how the probability distribution functions of pinning strengths β depend on properties of the network topology (Fig. 6c), we observe in ER networks that both the mean value and the variance of the distribution increase with increasing the wiring probability. In contrast, control of SF networks requires smaller pinning strengths.
As with the networks of second-order Kuramoto oscillators, we observe again that the pinning strength β is highly correlated to the node degree k, for both SF and ER networks (Fig. 6d,e), however, the correlation points here to a nonlinear relationship. The averaged values of degree k of driver nodes for both networks are 〈 〉  k 3, however, the averaged pinning strength β are β 〈 〉  260 and β 〈 〉  540 for ER and SF networks, respectively. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
Complex systems as diverse as the brain, power grids, food webs, social systems, and even the climate are often modeled as networks of dynamical units, and control of dynamical processes on such networks with complex topologies remains a highly active interdisciplinary field of research 4,[53][54][55] . Control strategies that are solely based on the network's topological properties 4 disregard the system's dynamics, which is usually highly nonlinear. Addressing this issue, we here proposed a straightforward topology-dynamics-based control approach for multi-dimensional complex networked dynamical systems (with and without self-dynamics). Our approach allows to find a suitable control function as well as a finite set of driver nodes (in polynomial time) to steer the system to the desired state. Moreover, the iterative scheme we developed to identify driver nodes can be combined with other control strategies such as the master stability function formalism (MSF) 26 , as discussed in detail in Supplementary Information. In contrast to MSF, our proposed approach can be used to determine the control function and the pinning strength(s) required for controlling complex systems that consist of non-identical components. This makes our approach different from most of the current control strategies.
By extensive numerical simulations and analyses of different multi-dimensional dynamics (second-order Kuramoto and Rössler oscillators) on different networks, we observed properties of driver nodes, such as their degree k and the pinning strength β to be related to both the underlying coupling topology and dynamics. For instance, our results indicate that in case of second-order Kuramoto dynamics the pinning density n D mainly depends on the network's mean degree while in the case of Rössler oscillators, it depends on both the network's mean degree and size.
Our approach relies on the existence of a fixed point of the system as a desired state, along with full observability and accessibility of the system. We have shown that one can employ our control approach also in cases where the desired state is a nominal trajectory by combining it with the master stability function formalism. It would be interesting to re-express our driver-node-identification scheme using the induced (matrix) norm 10,11,56 .
So far we have demonstrated our control approach using exemplary multi-dimensional dynamical systems, whose equations of motion are known, and using some prototypical network topologies. In a next step, our approach needs to be supplemented with a suitable minimisation of the control energy 55,57 derived from the corresponding dynamics. Once this has been achieved, we envision data-driven approaches to construct both dynamics and relevant aspects of the coupling 49,[58][59][60] . For these approaches, however, one should take into account that a fully observable system might have poor observability properties 61,62 . Nevertheless, such approachestogether with recently proposed data-driven methods to assess resilience 63 -would enable control of complex natural [64][65][66][67] and man-made systems 68,69 , where both topology and dynamics play essential roles in their collective behavior. Basic idea to determine the control function. Let us describe a general formalism which can be used to determine a control function for given coupled dynamics. Here, we assume that the system has two-dimensional state variables and that the dynamics is given by

Methods
with f and g being arbitrary and non-linear functions. We assume that the coupled dynamics has a fixed point at (x * , y * ) (by equating f(x, y) = g(x, y) = 0). The Jacobian matrix evaluated at the fixed point reads x x y y x x y y x x y y x x y y The center c 1 and radius r 1 of first the Gershgorin disk of this matrix read x x y y If r 1 + c 1 ≥ 0, this disk will have some part in the positive half plane of complex plane and consequently, it might cause a positive real part in eigenvalues. Here, we assume that this situation happens, and we determine a control function that allows stabilising the system's dynamics around its fixed point. We start with a transformation that enables steering the first Gershgorin disk to the left half plane, and the only disk which remains to be controlled will be the second one. We use . Accordingly, the dynamics of transformed variables read with the new Jacobian matrix evaluated at the fixed point, x x y y x x y y x x y y x x y y One can easily find elements of this matrix to check whether the first Gershgorin disk is placed in the left half plane or not. One finds where u is the second argument of f and g, which is in fact the variable y. Now, the new first Gershgorin disk is centered at  c 1 and has radius  r 1 x x y y x x y y It s hould be noted that adding the term β − −   ⁎ y y ( ) to Eq. 15 is equivalent to adding −β(y − y * ) − hβ(x − x * ) to Eq. 8.

Data availability
The data and code utilized in this study are available from the corresponding author upon reasonable request.