## Abstract

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.

## Introduction

Emergent behaviour of complex systems can be either desired (e.g. learning or intelligence) or undesired (e.g. hurricanes or cascading failures)^{1,2}. Generally, one aims to retain complex systems closer to the desired state, and control of such systems is important in order to steer the system’s dynamics to that state (assuming that it exists)^{3,4,5}. Contemporary control strategies can be divided into four major categories, depending on the control objective and the type of system dynamics they are suited for.

- (i)
*Linear systems or linearized dynamics control*aims to fully control linear systems or general non-linear systems near a nominal trajectory or a fixed point. Based on this objective, proposed strategies are based on Kalman rank condition^{6}, structural controllability^{7,8,9}, structural adaptation and contraction^{10,11}, exact controllability^{12}, path reachability^{13,14}, and on spectral properties of the Jacobian matrix^{15}. - (ii)
*Non-linear systems control*initially pursued the ambitious full controllability of general non-linear systems^{16,17,18,19}but was proven to be NP-hard for some cases^{20}, which led to weaker notions such as local accessibility/controllability^{20}. - (iii)
*Final-state-oriented control*does not actually intend to fully control the state of the system, but aims to steer the system towards a desired long-term behavior. Control through adjustments of the initial condition of the dynamics (placing the system in a new basin of attraction)^{21}, control through change of parameters and exploiting network regulatory attributes^{22}and controlling chaos^{23}are known works conducted in this category. Another strategy developed recently in this category, based on the concepts of “feedback vertex set” and “determining nodes”, is the so called “feedback vertex control”^{24,25}which is almost independent of the underlying dynamics of the network and mostly depends on the wiring of the system and knowledge of the target invariant sets. - (iv)
*Collective behavior control*intends to control the collective complex behavior of a system. Some well-known strategies in this category are based on the master stability formalism^{26,27}and on pinning control^{28,29}, motivated by publications on chaos control^{30}, and nonlinear distributed control protocols^{31}aim at identifying the required pinning term to be added to the dynamics of all nodes.

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 \({\mathscr{O}}(\sqrt{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 \({\mathscr{O}}({N}^{2}\,{\mathrm{ln}}^{2}(N))\). 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 \({\mathscr{O}}{\mathrm{(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 *i*th 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 \({\nu }_{i}(t)={\dot{\theta }}_{i}(t)\), which fully determine the system’s dynamics. If the phase *θ*_{i}(*t*) of oscillator *i* approaches the phase of the mean field (i.e., \(\psi =arctan(\Im ({\sum }_{j}{e}^{i{\theta }_{j}})\Re ({\sum }_{j}{e}^{i{\theta }_{j}}))\), 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 \(-{\beta }_{i}({x}_{i}-{x}_{i}^{\ast })\), where *β*_{i} is the pinning strength and *x*_{i} the state variable of node *i*. Here \({x}_{i}^{\ast }\) 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.

Now, if one defines \({\dot{\theta }}_{i}={\nu }_{i}\) and \({\dot{\nu }}_{i}\) as two coupled dynamical equations, the first *N* disks (out of the 2*N* Gershgorin disks for the corresponding Jacobian matrix evaluated at the fixed point – which are associated with \(\partial {\dot{\theta }}_{i}\) (here ∂ refers to partial differentiations of \(\dot{\theta }\) with respect to *ν* or *θ*) – will be centered at the origin and will have unit radius each (Supplementary Information). Eigenvalues of the Jacobian matrix, which lie on the union of all Gershgorin disks, thus may possess a positive real part. This can be avoided by a change of variables by the canonical transformations \({\tilde{\theta }}_{i}={\theta }_{i}\) and \({\tilde{\nu }}_{i}={\nu }_{i}+\kappa {\theta }_{i}\), with *κ* = 1, and by rewriting the set of Eq. (1). The new governing equations then read:

The Gershgorin disks related to \(\partial {\dot{\tilde{\nu }}}_{i}\) of the corresponding Jacobian (Supplementary Information) have the following centers *C*_{i} and corresponding radii *R*_{i}:

for *i* = *N* + 1, *N* + 2, …, 2*N*. Nodes for which the value of *C*_{i} + *R*_{i} is positive are “problematic” since their dynamics can induce a positive real part of eigenvalues. It can be shown (Methods) that adding a term of the form \(-{\beta }_{i}({\tilde{\nu }}_{i}-{\tilde{\nu }}_{i}^{\ast })\) to the r.h.s. of Eq. (2) for \({\dot{\tilde{\nu }}}_{i}\) or, equivalently, a term of the form \(-{\beta }_{i}({\nu }_{i}-{\nu }_{i}^{\ast })-{\beta }_{i}({\theta }_{i}-{\theta }_{i}^{\ast })\) to the equation for \({\dot{\nu }}_{i}\), enables us to move the eigenvalues of the Jacobian to the left half-plane of the complex plane. This possible control strategy can be interpreted as pinning the variables to the values of their associated components of the desired state, which assures the stability of such a state if the pinning strength *β*_{i} for node *i* satisfies the following condition:

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}.Step 2: repetition of step 1 separately for each of the remaining

*N*− 1 nodes leads to the set of values {*e*_{1},*e*_{2}, …,*e*_{N}}.

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 \(\simeq 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 (\(\beta \simeq 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 \(\beta \simeq 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 \(\langle \omega \rangle \simeq 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 (\(\langle k\rangle \simeq 4\)). These findings demonstrate that properties of driver nodes – namely their natural frequency and their mean degree – are related to the underlying coupling topology.

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 〈\(\langle k\rangle \simeq 11\) and \(\langle \beta \rangle \simeq 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 \(-{\beta }_{i}({\tilde{\nu }}_{i}-{\tilde{\nu }}_{i}^{\ast })\) 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

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

with

**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 \({g}_{ij}={\delta }_{ij}({\sum }_{l\mathrm{=1}}^{N}{a}_{il})-{a}_{ij}\), with *a*_{ij} denoting elements of **A**.

To control this network, we derive a term (Supplementary Information) of the form \(-{\beta }_{i}({x}_{3}^{i}-{x}_{3}^{i\ast })+{\beta }_{i}{\alpha }_{i}\)\(({x}_{1}^{i}-{x}_{1}^{i\ast })+{\beta }_{i}{\alpha }_{i}\mathrm{(1}+a)({x}_{2}^{i}-{x}_{2}^{i\ast })\) that – when added to the equation for \(\frac{d}{dt}{{x}_{3}}^{i}\) (see Eq. (6)) – enables us to move the eigenvalues of the corresponding Jacobian to the left half-plane of the complex plane. Here \({x}_{j}^{i\ast }\) (*j*∈{1, 2, 3}) are the stationary solutions of Eq. (6). Also, \({\alpha }_{i}\mathrm{=(2}+a\mathrm{)[2}+\sigma {k}_{i}^{{\rm{in}}}]\), and \({k}_{i}^{{\rm{in}}}={\sum }_{j}{a}_{ij}\) is the sum of weights connected to node *i*. The pinning strength *β*_{i} of node *i* satisfies the following relation

where the term \({\mathscr{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}^{{\rm{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.

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 \(\langle k\rangle \simeq 3\), however, the averaged pinning strength *β* are \(\langle \beta \rangle \simeq 260\) and \(\langle \beta \rangle \simeq 540\) for ER and SF networks, respectively.

## 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 approaches – together 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.

## Methods

### Linear stability, spectrum of eigenvalues, and Gershgorin circle theorem

A system is linearly stable around a fixed point, if all eigenvalues of the corresponding Jacobian matrix have a negative real part. The vanishing of the positive real part can be verified using the Gershgorin circle theorem^{39}, which states that all eigenvalues of a given *N* × *N* matrix M lie within the union of closed Gershgorin disks *D*_{i} with *i* = 1, …, *N*. The *i*th disk’s center is positioned at *C*_{i} = *M*_{i,i} and its radius is \({R}_{i}={\sum }_{j\ne i}|{M}_{i,j}|\), where *M*_{i,j} is an entry of M. Accordingly, to ensure that all eigenvalues lie in the negative half-plane of the complex plane (i.e., all components of dynamical system will approach the fixed point), all Gershgorin disks should be moved to the left half-plane of the complex plane.

### 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

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

The center *c*_{1} and radius *r*_{1} of first the Gershgorin disk of this matrix read

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

where \(h={\rm{sgn}}(\frac{\partial f}{\partial y}{|}_{\mathop{y={y}^{\ast }}\limits^{x={x}^{\ast }}})\). Accordingly, the dynamics of transformed variables read

with the new Jacobian matrix evaluated at the fixed point,

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 \({\tilde{c}}_{1}\) and has radius \({\tilde{r}}_{1}\)

It is obvious that \({\tilde{c}}_{1}+{\tilde{r}}_{1}=0\) holds if \({\tilde{c}}_{1}\) is negative which is true in practical cases.

Next, we should make sure that the second Gershgorin disk has no positive part. To do so, we first calculate its center \({\tilde{c}}_{2}\) and radius \({\tilde{r}}_{2}\)

Now, if \({\tilde{c}}_{2}+{\tilde{r}}_{2}\ge 0\), one can add a pinning term of the form of \(-\beta (\tilde{y}-{\tilde{y}}^{\ast })\) to the r.h.s. of Eq. 15. The pinning term adds a (−*β*) to the r.h.s. of Eq. 21 and moves the center to the left. To make sure that the Gershgorin disk has no positive real part, *β* should satisfy the f ollowing condition

It s hould be noted that adding the term \(-\beta (\tilde{y}-{\tilde{y}}^{\ast })\) to Eq. 15 is equivalent to adding −*β*(*y* − *y*^{*}) − *hβ*(*x* − *x*^{*}) to Eq. 8.

### Control function for networks of second-order kuramoto oscillators

The equation of motion of the network of second-order Kuramoto oscillators (Eq. (1)) can be rewritten as

for which the Jacobian matrix evaluated at the fixed point reads

*I* is the *N* × *N* unit matrix and entries of matrix M are

It is evident that the first *N* (out of 2*N*) Gershgorin dis ks for this Jacobian are centered at the origin and ha ve unit radii. The corresponding eigenvalues thus may possess a positive real part. To avoid this, we define tw o new variables, \({\tilde{\theta }}_{i}={\theta }_{i}\) and \({\tilde{\nu }}_{i}={\nu }_{i}+{\theta }_{i}\). This change of variables changes the g overn ing equations to (see Eq. (2)):

for which the Jacobian matrix evaluated at the fixed point has the following form

The entries of matrix M are

Here, \({\theta }_{i}^{\ast }\) denotes the equilibrium (or resting) point of node *i*, i.e., the value of the state variable for which \({\dot{\tilde{\theta }}}_{i}={\dot{\tilde{\nu }}}_{i}=0\), and it may be a stable or an unstable fixed point. It is easy to see that adding a term of the form \(-{\beta }_{i}({\tilde{\nu }}_{i}-{\tilde{\nu }}_{i}^{\ast })\) to Eq. (2) enables us to move the eigenvalues of the Jacobian to t he left half-plane of the complex plane, if the pinning strength *β*_{i} satisfies condition (4).

## Data availability

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

## References

- 1.
Haken, H.

*Synergetics - An Introduction and Advanced Topics*(Springer, Berlin, 2004). - 2.
Kwapień, J. & Drożdż, S. Physical approach to complex systems.

*Phys. Rep.***515**, 115–226 (2012). - 3.
Motter, A. E. Networkcontrology.

*Chaos***25**, 097621 (2015). - 4.
Liu, Y.-Y. & Barabási, A.-L. Control principles of complex systems.

*Rev. Mod. Phys.***88**, 035006 (2016). - 5.
Barabási, A. & Posfai, M.

*Network Science*, 1st edn (Cambridge University Press, Cambridge, UK, 2016). - 6.
Kalman, R. E. Mathematical description of linear dynamical systems.

*J. SIAM Series A Control***1**, 152–192 (1963). - 7.
Lin, C.-T. Structural controllability.

*IEEE Trans. Autom. Control***19**, 201–208 (1974). - 8.
Liu, Y.-Y., Slotine, J.-J. & Barabási, A.-L. Controllability of complex networks.

*Nature***473**, 167–173 (2011). - 9.
Whalen, A. J., Brennan, S. N., Sauer, T. D. & Schiff, S. J. Observability and controllability of nonlinear networks: The role of symmetry.

*Phys. Rev. X***5**, 011005 (2015). - 10.
Lohmiller, W. & Slotine, J.-J. E. On contraction analysis for non-linear systems.

*Automatica***34**, 683–696 (1998). - 11.
DeLellis, P., Di Bernardo, M., Gorochowski, T. E. & Russo, G. Synchronization and control of complex networks via contraction, adaptation and evolution.

*IEEE Circ. Syst. Mag.***10**, 64–82 (2010). - 12.
Yuan, Z., Zhao, C., Di, Z., Wang, W.-X. & Lai, Y.-C. Exact controllability of complex networks.

*Nat. Commun.***4**, 2447 (2013). - 13.
Parlangeli, G. & Notarstefano, G. On the reachability and observability of path and cycle graphs.

*IEEE Trans. Autom. Control***57**, 743–748 (2012). - 14.
Bai, Y.-N., Wang, L., Chen, M. Z. Q. & Huang, N. Controllability emerging from conditional path reachability in complex networks.

*Int. J. Robust. Nonlinear Control*(2017). - 15.
Skardal, P. S. & Arenas, A. Control of coupled oscillator networks with application to microgrid technologies.

*Sci. Adv.***1**, e1500339 (2015). - 16.
Elliott, D. L. A consequence of controllability.

*J. Differ. Equ.***10**, 364–370 (1971). - 17.
Rugh, W. J.

*Linear System Theory, Information and System Sciences*(Prentice Hall, Englewood Cliffs, NJ, USA, 1993). - 18.
Haynes, G. W. & Hermes, H. Nonlinear controllability via Lie theory.

*SIAM J. Control***8**, 450–460 (1970). - 19.
Sontag, E. D.

*Mathematical control theory: deterministic finite dimensional systems*, vol. 6 (Springer Science & Business Media, 2013). - 20.
Sontag, E. D. Controllability is harder to decide than accessibility.

*SIAM J. Control Optim.***26**, 1106–1118 (1988). - 21.
Cornelius, S. P., Kath, W. L. & Motter, A. E. Realistic control of network dynamics.

*Nat. Commun.***4**, 1942 (2013). - 22.
Wang, L.-Z.

*et al*. A geometrical approach to control and controllability of nonlinear dynamical networks.*Nat. Commun.***7**, 11323 (2016). - 23.
Ott, E., Grebogi, C. & Yorke, J. A. Controlling chaos.

*Phys. Rev. Lett.***64**, 1196–1199 (1990). - 24.
Mochizuki, A., Fiedler, B., Kurosawa, G. & Saito, D. Dynamics and control at feedback vertex sets. II: A faithful monitor to determine the diversity of molecular activities in regulatory networks.

*J. Theor. Biol.***335**, 130–146 (2013). - 25.
Za˜nudo, J. G. T., Yang, G. & Albert, R. Structure-based control of complex networks with nonlinear dynamics.

*Proc. Natl. Acad. Sci. (USA)***114**, 7234–7239 (2017). - 26.
Pecora, L. M. & Carroll, T. L. Master stability functions for synchronized coupled systems.

*Phys. Rev. Lett.***80**, 2109–2112 (1998). - 27.
Sorrentino, F., di Bernardo, M., Garofalo, F. & Chen, G. Controllability of complex networks via pinning.

*Phys. Rev. E***75**, 046103 (2007). - 28.
Li, X., Wang, X. & Chen, G. Pinning a complex dynamical network to its equilibrium.

*IEEE Trans. Circ. Syst. I: Regular Papers***51**, 2074–2087 (2004). - 29.
Yu, W., Chen, G., Lu, J. & Kurths, J. Synchronization via pinning control on general complex networks.

*SIAM J. Control Optim.***51**, 1395–1416 (2013). - 30.
Boccaletti, S., Grebogi, C., Lai, Y.-C., Mancini, H. & Maza, D. The control of chaos: theory and applications.

*Phys. Rep.***329**, 103–197 (2000). - 31.
Monteil, J. & Russo, G. On the design of nonlinear distributed control protocols for platooning systems.

*IEEE Control Syst. Lett.***1**, 140–145 (2017). - 32.
Pequito, S., Preciado, V. M., Barabási, A.-L. & Pappas, G. J. Trade-offs between driving nodes and time-to-control in complex networks.

*Sci. Rep.***7**, 39978 (2017). - 33.
Hopcroft, J. E. & Karp, R. M. An

*n*5=2 algorithm for maximum matchings in bipartite graphs.*SIAM J. Comput.***2**, 225–231 (1973). - 34.
Hautus, M. L. J. Controllability and observability conditions of linear autonomous systems.

*Proceedings of the Koninklijke Nederlandse Akademie Van Wetenschappen Series a-Mathematical Sciences***72**, 443 (1969). - 35.
Nacher, J. C. & Akutsu, T. Dominating scale-free networks with variable scaling exponent: heterogeneous networks are not difficult to control.

*New J. Physics***14**, 073005 (2012). - 36.
Nacher, J. C. & Akutsu, T. Structural controllability of unidirectional bipartite networks.

*Sci. Rep***3**, 1647 (2013). - 37.
Fomin, F. V., Grandoni, F. & Kratsch, D. A measure & conquer approach for the analysis of exact algorithms.

*J. ACM***56**, 25 (2009). - 38.
Arenas, A., Dìaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou, C. Synchronization in complex networks.

*Phys. Rep.***469**, 93–153 (2008). - 39.
Golub, G. H. & Van Loan, C. F.

*Matrix computations*, third edn (John Hopkins University Press, Baltimore, MD, USA, 1996). - 40.
Erdős, P. & Rényi, A. On the evolution of random graphs.

*Publ. Math. Inst. Hung. Acad. Sci.***5**, 17–61 (1960). - 41.
Albert, R. & Barabási, A.-L. Statistical mechanics of complex networks.

*Rev. Mod. Phys.***74**, 47–97 (2002). - 42.
Grigg, C.

*et al*. The IEEE reliability test system-1996. A report prepared by the reliability test system task force of the application of probability methods subcommittee.*IEEE Trans. Power Syst.***14**, 1010–1020 (1999). - 43.
Dörfler, F., Chertkov, M. & Bullo, F. Synchronization in complex oscillator networks and smart grids.

*Proc. Natl Acad. Sci. (USA)***110**, 2005–2010 (2013). - 44.
Apt, J. The spectrum of power from wind turbines.

*J. Power Sources***169**, 369–374 (2007). - 45.
Baile, R. & Muzy, J.-F. Spatial intermittency of surface layer wind fluctuations at mesoscale range.

*Phys. Rev. Lett.***105**, 254501 (2010). - 46.
Wood, R. & Field, P. R. The distribution of cloud horizontal sizes.

*J. Climate***24**, 4800–4816 (2011). - 47.
Milan, P., Wächter, M. & Peinke, J. Turbulent character of wind energy.

*Phys. Rev. Lett.***110**, 138701 (2013). - 48.
Tabar, M. R. R.

*et al*. Kolmogorov spectrum of renewable wind and solar power fluctuations.*Eur. Phys. J. Spec. Top.***223**, 2637–2644 (2014). - 49.
Anvari, M., Tabar, M. R. R., Peinke, J. & Lehnertz, K. Disentangling the stochastic behavior of complex time series.

*Sci. Rep.***6**, 35435 (2016). - 50.
Gambuzza, L. V., Buscarino, A., Fortuna, L., Porfiri, M. & Frasca, M. Analysis of dynamical robustness to noise in power grids.

*IEEE Trans. Emerg. Sel. Topics Circuits Syst.***7**, 413–421 (2017). - 51.
Schäfer, B.

*et al*. Escape routes, weak links, and desynchronization in fluctuation-driven networks.*Phys. Rev. E***95**, 060203 (2017). - 52.
Lehnertz, K., Zabawa, L. & Tabar, M. R. R. Characterizing abrupt transitions in stochastic dynamics.

*New J. Physics***20**, 113043 (2018). - 53.
Gates, A. J. & Rocha, L. M. Control of complex networks requires both structure and dynamics.

*Sci. Rep.***6**, 24456 (2016). - 54.
Kim, J. Z.

*et al*. Role of graph architecture in controlling dynamical networks with applications to neural systems.*Nat. Phys.***14**, 91 (2018). - 55.
Lindmark, G. & Altafini, C. Minimum energy control for complex networks.

*Sci. Rep.***8**, 3188 (2018). - 56.
Vidyasagar, M.

*Nonlinear systems analysis*(SIAM, Philadelphia, 2002). - 57.
Yan, G., Ren, J., Lai, Y.-C., Lai, C.-H. & Li, B. Controlling complex networks: How much energy is needed?

*Phys. Rev. Lett.***108**, 218703 (2012). - 58.
Rahimi Tabar, M. R.

*Analysis and Data-Based Reconstruction of Complex Nonlinear Dynamical Systems: Using the Methods of Stochastic Processes*(Springer, Cham-Switzerland, 2019). - 59.
Prusseit, J. & Lehnertz, K. Measuring interdependences in dissipative dynamical systems with estimated Fokker-Planck coefficients.

*Phys. Rev. E***77**, 041914 (2008). - 60.
Friedrich, R., Peinke, J., Sahimi, M. & Tabar, M. R. R. Approaching complexity by stochastic methods: From biological systems to turbulence.

*Phys. Rep.***506**, 87–162 (2011). - 61.
Aguirre, L. A. & Letellier, C. Observability of multivariate differential embeddings.

*J. Phys. A: Mathematical and General***38**, 6311–6326 (2005). - 62.
Aguirre, L. A., Portes, L. L. & Letellier, C. Structural, dynamical and symbolic observability: From dynamical systems to networks.

*Plos One***13**, 1–21 (2018). - 63.
Rings, T.

*et al*. Traceability and dynamical resistance of precursor of extreme events.*Sci. Rep.***9**, 1744 (2019). - 64.
Li, A., Cornelius, S. P., Liu, Y.-Y., Wang, L. & Barabási, A.-L. The fundamental advantages of temporal networks.

*Science***358**, 1042–1046 (2017). - 65.
Tang, E. & Bassett, D. S. Colloquium: Control of dynamics in brain networks.

*Rev. Mod. Phys.***90**, 031003 (2018). - 66.
Tu, C.

*et al*. Warnings and caveats in brain controllability.*NeuroImage***176**, 83–91 (2018). - 67.
Matamalas, J. T., Arenas, A. & Gòmez, S. Effective approach to epidemic containment using link equations in complex networks.

*Sci. Adv*. 4 (2018). - 68.
Pagani, G. A. & Aiello, M. The power grid as a complex network: a survey.

*Physica A: Statistical Mechanics and its Applications***392**, 2688–2700 (2013). - 69.
Tang, Y., Gao, H., Zhang, W. & Kurths, J. Leader-following consensus of a class of stochastic delayed multi-agent systems with partial mixed impulses.

*Automatica***53**, 346–354 (2015).

## Acknowledgements

The authors would like to thank Thorsten Rings for interesting discussions and for critical comments on earlier versions of the manuscript. This work was supported by the Verein zur Foerderung der Epilepsieforschung e.V. (Bonn). MRRT acknowledges partial support by INSF under Grant No. 96000289.

## Author information

### Affiliations

### Contributions

M.B., M.R.R.T., J.P. and K.L. designed the study. K.L., M.R.R.T. and M.B. wrote the manuscript. Computer codes were written by M.B. (networks of Kuramoto oscillators, IEEE power grid), H.A. (networks of Rössler oscillators, master stability formalism, feedback vertex set), and T.M. All authors reviewed the manuscript.

### Corresponding author

## Ethics declarations

### Competing interests

The authors declare no competing interests.

## Additional information

**Publisher’s note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

## Rights and permissions

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

## About this article

### Cite this article

Bahadorian, M., Alimohammadi, H., Mozaffari, T. *et al.* A topology-dynamics-based control strategy for multi-dimensional complex networked dynamical systems.
*Sci Rep* **9, **19831 (2019). https://doi.org/10.1038/s41598-019-56259-4

Received:

Accepted:

Published:

## Comments

By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.