The complex networks we envision controlling are the common case of nonlinear, dissipative dynamical systems, where the time evolution of the component dynamical variables is governed by ordinary differential equations (ODEs). The roles of the individual nodes and the structure of the network as a whole are assumed to be reflected in these ODEs. For example, if there is a link in the network from node i to node j, one expects a term in the j-th differential equation involving the i-th state variable.Thus for the operational purpose of “control”, we can completely characterize a given complex network by:
• a (vector) function “f” that gives the “right-hand side” of the first-order ODEs defining the dynamics;
• two (vector) functions “g” and “h” that define, respectively, the inequality and equality constraints on the eligible interventions.
Note that from this point onward, we are no longer talking about NECO in the abstract, but rather in regard to its software implementation in Python. Thus f, g, and h here denote Python functions that are to be implemented by the user and input to neco. These should be regarded as distinct from the corresponding mathematical abstractions in Fig. 1, which are similarly labeled as F, g, and h (in Roman font) to facilitate comparison. In particular, the Python functions will generally take different arguments, and in a different order, than their abstract counterparts. For example, f—the Python implementation of the derivatives function F—can accept a number of optional parameters that are needed to compute the right-hand side of the ODEs, parameters that were omitted for notational clarity in Fig. 1.
Also, note that within the source code and from this point onward in this protocol, we use the letter 'y' instead of 'x' to denote the state of a dynamical system, to follow the library SciPy's notational conventions for systems of ODEs.
Defining the ODE system:
To apply this software to find compensatory perturbations in a given dynamical system, one needs only to supply a Python function of the form
f(y, t, param1, param2, …),
which, given a system state y and optional parameters param1, param2, etc., evaluates the right-hand side of the ordinary differential equations that govern the time evolution of the system. We anticipate that the user will frequently want to use this function for other purposes, e.g., integration of system orbits for the purposes of visualization. Thus, we have intentionally required the call signature of f to take the form above so that it can be used with SciPy's tools for ODEs (in particular their numerical integrator, odeint) without modification. Note, however, that neco assumes that the ODE system described by f be autonomous—that is, it has no explicit time dependence. Hence, even though the time variable t appears as a required second argument, it should not actually be used within the calculation of the return value of f. The condition of being autonomous does not sacrifice any generality, since every non-autonomous system of ODEs (those with an explicit dependence on t) can be transformed into an equivalent autonomous system by augmenting the state space with an additional “dummy” variable representing t, whose time derivative is then the constant 1.
Note that we impose no specific conditions on the nature of the dynamics defined by f; notably, it can contain arbitrary nonlinearities and be arbitrarily high-dimensional. The only (mild) requirement is that f be differentiable with respect to the system state y, so that a Jacobian matrix can be calculated.
Defining constraints on the eligible interventions:
Given the network state y0 at a particular time, one can specify the subset of network states y that are reachable through eligible perturbations from y0 using vector constraint functions of the form g(y) and h(y) for inequality and equality constraints, respectively. Each of these functions operates by taking a state y of the network and returning a NumPy array with lengths equal to the number of inequality or equality constraints, respectively. States that can be reached by eligible control perturbations are then those that satisfy g(y) <= 0 and h(y) == 0, where the inequality and equality are taken to apply individually to every element of the returned arrays.
A common occurrence of inequality/equality constraints is one in which some dynamical variables (components of y) are allowed to be perturbed over a certain range and/or cannot be perturbed at all. This special case of component-wise constraints can always be expressed using g and h above, but, for convenience, we have provided the option to specify vectors of lower (lb) and upper (ub) bounds on eligible perturbed states y, which must then satisfy
lb[i] <= y[i] <= ub[i],
for every component i. This includes in particular the case in which a state variable i cannot be perturbed at all, i.e., lb[i] = ub[i] = y0[i].
Each of lb, ub, g, and h is an optional keyword argument for the function neco. If either vector lb or vector ub is not supplied, the coordinates of the eligible states are assumed to be unbounded from below or above, respectively, before the imposition of any constraints specified by g(y) and/or h(y). Individual components can be made unbounded with respect to lb and ub by filling the appropriate positions in these vectors with the NumPy constants -inf and +inf, respectively.
Parameters for the control procedure:
The enclosed implementation of NECO accepts a handful of parameters that govern its internal operation. For the motivation behind each of these parameters, as well as guidance in how to choose their values based on the time and length scales of a particular dynamical system, we refer the reader to ref. . Each of these parameters is implemented as a keyword argument in the neco Python function, which we summarize here:
eps0: Minimum size of the incremental perturbation at every iteration (eps0 > 0), to ensure that the algorithm makes a non-zero step every time.
eps1: Maximum size of the incremental perturbation at every iteration (eps1 > eps0), to ensure that the forecasted perturbation using the variational matrix is valid (at both the initial time and the forecasted time of closest approach to the target).
it_max: Maximum number of iterations of the algorithm. If this number is exceeded before finding a compensatory perturbation, the function will return 1, indicating failure.
t_max: Time window over which to search for the closest approach to the target at every iteration.
dt: Integration time step.
t_test: Time window over which to test convergence of the system's orbit to the target state (generally t_test >> t_max)
tol: Numerical tolerance for convergence to the target state.
Thus, if a state y evolves to within a ball of radius tol around the target state, and within t_test time units, it is considered to be in the basin of attraction of the target and the algorithm declares success.
Other optional parameters:
jac(y, t, param1, param2, …): Function that computes the Jacobian matrix J of f as a 2D NumPy array, where J[i, j] = df[i]/dx[j]. If not supplied, the Jacobian will be calculated from f numerically.
dist(y1, y2): Function that defines the distance between two states, y1 and y2, in the system under consideration. By default, the Euclidean norm is used.
n_test: Every n_test iterations, test whether the current initial state attracts to the target. By default, n_test = 1.
full_output: If full_output is False, the method simply returns the final perturbed state. If full_output is True, the method returns a dictionary containing more detailed information (see "Return values" below). This option is False by default.
args: A tuple of additional parameters (other than y and t) that will be passed to the derivatives function f and the Jacobian function jac (if supplied).
(if full_output is False)
y0_prime: The final perturbed state.
(if full_output is True)
info: Dictionary containing the following:
y0: Sequence of intermediate perturbed states at every iteration. The last element, y0[-1], is then the final perturbed state.
status: The value 0 indicates success, meaning that y0[-1] is an eligible state in the target's
basin of attraction. The value 1 indicates that the maximum iteration limit was exceeded.
n_iter: Total number of iterations taken before completion.
time: Total run time of the method (in seconds).
t_int, t_var, t_opt: Lists of run times (in seconds) taken at each iteration for three substeps
(integrating the system equations, integrating the variational equation, and forecasting the optimal incremental perturbation, respectively).
How to invoke neco:
Place the file named neco.py in a directory contained in your PYTHONPATH environment
variable. From that point, you will be able to import the control routine from the module of the same name:
from neco import neco
…within any Python program. Once the differential equations defining the desired network have
been implemented as a Python function f (as described in Procedure above), the control method can then be invoked within your script according to
output = neco(y0, yt, f, ****kwargs)
…where y0 and yt are the initial and desired target states of the network, respectively, and ****kwargs represents any (optional) additional arguments that one wishes to specify or change from their default values (constraints, the Jacobian matrix, iteration limits, and so on as described in the Procedure section above). For example, to invoke neco with minimum and maximum incremental perturbation sizes of 10-3 and 10-2, respectively, and an integration time window of 10 time units, one would enter
output = neco(y0, yt, f, eps0=0.001, eps1=0.01, t_max=10)
...and so on.
As part of this protocol we include the source code (necoexample.py) for two examples of finding control perturbations subject to constraints in a two-dimensional system defined by a particle under the influence of a 1D potential plus dissipation. These control problems are illustrated in Fig. 2 of ref. 1. These examples are provided mainly for the source code itself (in particular, to give guidance to those unfamiliar with Python on how to define a system of ODEs and call _neco with various parameters). To actually run the examples, simply invoke python on the example file from the command line