A machine learning based control of complex systems

TL;DR

This paper presents a novel machine learning-based control method for complex systems, utilizing symbolic dynamics and homoclinic orbits to steer chaotic systems towards desired states, demonstrated through simulations.

Contribution

It introduces a data-driven control strategy that does not rely on explicit system modeling but uses machine learning to approximate control signals based on current states.

Findings

Effective control of chaotic systems demonstrated in simulations.

Control strategy applicable to discrete maps, differential equations, and coupled networks.

Shows potential of machine learning in complex system regulation.

Abstract

In this work, inspired in the symbolic dynamic of chaotic systems and using machine learning techniques, a control strategy for complex systems is designed. Unlike the usual methodologies based on modeling, where the control signal is obtained from an approximation of the dynamic rule, here the strategy rest upon an approach of a function, that from the current state of the system, give the necessary perturbation to bring the system closer to a homoclinic orbit that naturally goes to the target. The proposed methodology is data-driven or can be developed in a based-model context and is illustrated with computer simulations of chaotic systems given by discrete maps, ordinary differential equations and coupled maps networks. Results shows the usefulness of the design of control techniques based on machine learning and numerical approach of homoclinic orbits.

Equations

$$s_{t+1}=f_r(s_t)+u_t,$$

$$f_{r}^{-1}(s_{t},\Sigma)=\left\{\begin{array}[]{ccc}f_{r,0}^{-1}(s_{t})&if&\sigma_{t}=0\\ f_{r,1}^{-1}(s_{t})&if&\sigma_{t}=1\\ \end{array}\right.,$$

$$s_0=f_{r,{\sigma }_1}^{-1}\circ \cdots \circ f_{r,{\sigma }_{l-1}}^{-1}\circ f_{r,{\sigma }_l}^{-1}(s^{\ast }+\delta )$$

$${\Sigma }_i={\sigma }_1^i,{\sigma }_2^i,\cdots ,{\sigma }_k^i,\underset{(l-k)times}{\underbrace{{\sigma }^{\ast },\cdots ,{\sigma }^{\ast }}}$$

$$s_0^i=f_r^{-l}(s_{R_i},{\Sigma }_i)$$

$$u_t=P(S_0,s_t)-s_t,$$

$$P(S_0,s_t)=argmi{n}_{s_0^i}\{|s_0^i-s_t|\}.$$

$$u\sim GP(m(s),k(s,s^{\prime })),$$

$$f_{r}^{-1}(s_{t},\Sigma)=\left\{\begin{array}[]{cc}\frac{1}{2}+\frac{\sqrt{1-s_{t}}}{2}&if~{}\sigma_{t}=0\\ \frac{1}{2}-\frac{\sqrt{1-s_{t}}}{2}&if~{}\sigma_{t}=1\\ \end{array}\right.,$$

$$x_{t+1}$$

$$y_{t+1}$$

$$\dot{x}(t)$$

$$\dot{y}(t)$$

$$\dot{z}(t)$$

$$s_{t+1}^i=f(s_t^i)-\frac{ϵ}{k_i}\sum _{j=1}^M{L}_{ij}f(s_t^j),$$

$${\mathbf{s}}_{t+1}=\mathbf{F}({\mathbf{s}}_t)=\mathbf{f}({\mathbf{s}}_t)-ϵ\mathbf{C}\mathbf{f}({\mathbf{s}}_t),$$

$${\bf A}=\left(\begin{array}[]{cccccc}0.9&0.&0.0333&0.0333&0.0333&0.\\ 0.&0.9&0.025&0.025&0.025&0.025\\ 0.0333&0.0333&0.9&0.&0.0333&0.\\ 0.0333&0.0333&0.&0.9&0.0333&0.\\ 0.02&0.02&0.02&0.02&0.9&0.02\\ 0.&0.05&0.&0.&0.05&0.9\\ \end{array}\right).$$

$$s_{t+1}^i=f(s_t^i)-\frac{ϵ}{k_i}\sum _{j=1}^M{L}_{ij}f(s_t^j)+u_t^i,$$

Full text

1 Introduction

There are experimental and/or numerical evidence of the chaotic behavior in systems from: population dynamics[11], epidemiology[31], ecosystems[29], cardiac rhythm[34], neurological systems[15], optoelectronic systems[14], chemical reactions[33], economic systems[2], communication systems[20] and a lot of more situations in physics, chemistry, biology and engineering[12]. Added to this, that behavior can be useful in areas such as secure communication, where chaotic systems offer an alternative for the development of new secure communication technologies [25, 30].

As can be seen in most of the before mentioned cases, it is desirable to have strategies that allow the evolution of the system to be regulated: the population of animal communities can be controlled to prevent extinctions or the spread of diseases, for obvious reasons the epidemic diseases must be controlled, the heart rate must be controlled in the case of arrhythmia, brain activity must be controlled in case of epilepsy or finally some chemical, physical or computational systems can be controlled for the benefit of all. This represents sufficient motivation to devote effort to the design models to forecast or control schemes for such systems.

The origin of chaotic behavior in dynamic systems can be explained, at least partially, from two different perspectives: the first is quantitative or statistical and is associated to the determination of quantities such as Lyapunov exponents that discriminate between chaotic or non-chaotic behavior. The second approach involves a geometric view for the dynamical system and attempts to isolate the geometrical structures that support the chaotic motion. The basic geometric objects, in this case, are the stable and unstable manifolds of fixed points of maps, or critical points or periodic orbits of ordinary differential equations. The homoclinic orbits are orbits that lies in both the stable and unstable manifolds of fixed points, or periodic orbits and particularly joins a saddle equilibrium point to itself.

The effect of the intersection of the stable and unstable manifolds, i. e., the homoclinic orbit, is nevertheless felt globally and offers a way in which simple local information can be extrapolated to complicated global behavior [16].

From a general perspective a dynamical system can be modeled as $s_{t}=f_{r}^{t}(s_{0})$, where $r$ is a set of parameters, $t\in\Re$ or $t\in\cal N$ represents the time, $s\in S$ the state of the system and $f$, a function that regulates the evolution of the system; so that, the control over them can only be achieved applying disturbances of some of these components.

When systems have chaotic behavior: on the one hand they have sensitivity to the initial conditions, which makes them difficult to predict in the long term, but on the other it offers the opportunity to change their behavior radically using small perturbations, this is, with a low cost and without producing important alterations in the system.

The seminal article on this topic presents a method for the control of chaotic systems, known as OGY method[21] and where the term *controlling chaos * was coined, is the archetype of the methods that belong to the class of those that disturb the parameter. As evidence of its importance, from the point of view of the basic research or the technology of the finding of Ott, Grebogi and York, there is the great amount of methods to stabilize chaotic systems that have been designed later, see example [10] and references therein. Another method almost as celebrated, as the previous one, is a feedback method proposed by Pyragas et. al. [23] belonging to the methods that disturb $f$. In this method, the stabilization of unstable periodic orbits of a chaotic system is achieved by combined feedback with the use of a specially designed external oscillator. Neither of the two methods require an a priory analytical knowledge of the dynamics system and may be applicable in experimental situations.

Among the control strategies of chaotic systems, there are a significant amount of control methods based on linear approximations of the dynamics around the target. When is it so, it is necessary that the linear approach is a good representation of the system or equivalently that the trajectory of the system is close enough to the target, which makes the control less efficient.

In this work, unlike those that model the dynamics in linear or nonlinear form, to design the control function, here a nonlinear control function is directly obtained by using Kernel Regression methods[26] and data from the observation of the system. The training set for the kernel machine, in this case, is constructed using a set of states that are close to points of a homoclinic orbit.

The article is organized as follow: in Section 2 a symbolic dynamic inspired but machine learning based approach to control of chaotic dynamic is presented. In Section 3, two versions of the solution to the problem are presented, depending on whether or not the symbolic dynamics of the system are known. Section 4 is devoted to give out some examples of the performance of the methodology and finally in Section 5 we give some concluding remarks.

2 A machine learning approach to control chaos

Colloquially, a control problem can be paraphrased as a sequence of three steps: (i) the identification of the control’s object, (ii) the selection of the control’s target and (iii) the design of a control strategy.

2.1 Control object and target

We will consider chaotic systems in one or more dimensions, continuous or discrete and represented by discrete maps, ordinary differential equations or coupled map networks. However, in order to show the ideas behind our control problem in a simple form, we will start considering as control object a chaotic, one-dimensional and uni-modal map $s_{t+1}=f_{r}(s_{t})$, for later in the results section use as control objects the rest of previously listed systems.

The targets, in all examples, are given by a unstable fixed-point solutions (${\bf s}^{*}$) of the dynamical system, although the control strategy can be adapted to the case of other types of periodic orbits.

2.2 Control strategy

Our approach to the control strategy propose to the estimate a small perturbation $u_{t}$, to the actual state $s_{t}$, in such a way that in the $(t+1)$-th iteration of $f$ the orbit come close to some point $s_{0}$, whose evolution has the form $s_{0},s_{1},\cdots,s_{k-1},s^{*}+\delta$, where $\delta$ is a small value. This is, we should to estimate a perturbation $u_{t}$, such that, if

[TABLE]

then $|s_{t+1}-s_{0}|<|s_{t}-s_{0}|$, with $s^{*}+\delta=f^{k}_{r}(s_{0})$.

In some sense, the idea behind this methodology is like to the OGY method, it consists in bringing the system closer, not immediately, to the control’s target but to points of orbits that naturally takes it as close as possible to target. However, unlike the OGY method, we do not disturb the parameters of the system or use explicitly a linear approximation of the system.

Thus, if we are able to estimate these new points and include it in the target, we expect that the control scheme to be more effective. The points that fulfill that condition, in the case of the fixed points of chaotic systems, are the points of homoclinic orbits. This orbits joins a saddle equilibrium point to itself[7], and offer a number of potential targets equal to the points in this orbit.

So, our control strategy start by determining the set of estates $\{s^{i}_{0}\}_{i=1}^{q}$. Clearly, these states can be estimated if we iterate the inverse function of $f$ (evolving in time reversal the system), from $s^{*}+\delta$, this is, $s_{0}=f^{-k}_{r}(s^{*}+\delta)$, but this can no be done simply because the inverse of a nonlinear function is multi-valued.

However, it is well known[4] that the phase space evolution of an uni-modal discrete map can be translated, in bi univocal way, into a binary representation by using a partition of the state space in two regions, labeled each of them by one symbol ($\sigma_{t}$), [math] or $1$. The symbolic representation is obtained replacing $s_{t}$ by the label associated to the element of the partition which belongs $s_{t}\to\sigma_{t}$. With this new information $f_{r}^{-1}$, can be written as as:

[TABLE]

where $f_{r,0}$ and $f_{r,1}$ are the monotone branch at the left or right from the maximum of $f$, respectively. In this form, every state $s_{t}$, can be represented by a sequence of symbols $\Sigma=\{\sigma_{j}\}_{j=1}^{l}$, given by the symbolic chain generated from the initial condition $s_{0}$. Here $l$ give the accuracy with which $s_{0}$ is represented. In particular, the states from which the system reaches a neighborhood of the fixed point are given by:

[TABLE]

For $l$ large enough, the iteration of (3) converges approximately to the real number $s_{0}$ independently of the value of $\delta$, if the sequence of symbols $\sigma_{t}$ is appropriate, i. e. if these sequences are of the form

[TABLE]

are large enough and $\sigma^{*}$ is the symbol associated to the interval that $s^{*}$ belongs.

The sequences $\Sigma_{i}$ form a set of $2^{k}$ different binary chains. Not all of them are admissible (see for example [4], Sec. 2.5.5), but the existence of some of them is guaranteed by the presence of homoclinic orbits in the attractor of $f$. This sequences can be used as a symbolic representation of segments of homoclinic orbits and will be used to estimate the $2^{k}$ initial conditions,

[TABLE]

with $s_{R_{i}}$ a random real number in the adequate domain, which places the system in a orbit towards the target $s^{*}$.

Now, in order to estimate the $u_{t}$ we start generating, from $q$ strings like (4), $q$ initial conditions $S_{0}=\{s^{1}_{0},s^{2}_{0},\cdots,s^{q}_{0}\}$, from which the system reach a neighborhood of $s^{*}$ in $l$ iterations.

With the before estimate states, the $N$ pairs $\{u_{t},s_{t}\}_{t=1}^{N}$ are generated using,

[TABLE]

from the elements of an typical orbit $\{s_{t}\}_{t=1}^{N}$ of the system. Here, $P(S_{0},s_{t})$ is a function that pick out, the element of $S_{0}$ nearest to $s_{t}$ is defined as:

[TABLE]

These pairs will later serve as training data for a nonlinear regression scheme $u_{t}=u(s_{t})$. In this way the control problem is turned into a regression problem, that we propose to solve using machine learning techniques.

As one, among many alternatives to approximate $u$, in this work we use a powerful set of strategies from machine learning known as Kernel Regression Methods [26, 13].

The main idea behind of these techniques is to map the observed examples (training set) into a new representation space (known as feature space), where object properties can be extracted more easily, using linear algorithms.

With the aim of showing how one of these techniques can be selected according to the particularities of the problem, in this case we present results with Gaussian Process[24] in the case of problems associated with a small number of data and Random Forest[8] in the case that the problem involve a large amount of data.

2.2.1 Gaussian Process Regression

Gaussian processes is a kernel-based Bayesian tool for modeling problems. Like other kernel-based methods, it combine a high flexibility of the model by working in often infinite-dimensional feature spaces with the simplicity that all operations are performed in the lower-dimensional input space using positive definite kernels.

Here we will understand these processes according to their usual definition: a set of indexed random variables, such that every finite subset of those random variables has a multivariate normal distribution. The distribution of a $GP$ is the joint distribution of all those random variables, and as such, it is a distribution over functions with a continuous domain.

These process are specified by its mean function $m(s)$ and covariance function $k(s,s_{i})$ and are a natural generalization of the Gaussian distribution whose mean and covariance are represented by a vector and a matrix, respectively. We will represent this process, with the usual notation, as:

[TABLE]

and we will read this, as: *the function $u$ is distributed as a Gaussian process with prior mean function $m(x)$ and covariance function $k(s,s^{\prime})$ *. If we generate a random vector from this distribution, $u\sim{\cal N}(m(s)$, $k(s,s^{\prime}))$, this vector will have the function values $u(s)$ as coordinates indexed by $s$.

This Gaussian process directly captures the model uncertainty but are computationally expensive in the case of great amount of data.

2.2.2 Random Forest Regression

Decision tree learning is one of the predictive modeling approaches used in statistics, data mining and machine learning. It uses a decision tree, as a predictive model, to go from observations about an item (represented in the branches) to response about the item’s target value (represented in the leaves).

Random forests or random decision forests are an ensemble randomized decision trees designed for classification or regression. This strategy can be seen as a kernel method [27, 6] that operate by constructing a multitude of decision trees at training time and average prediction of the individual trees.

This methodology runs efficiently on large databases and reduces overfitting in decision trees, helping us to improve the accuracy.

3 Model-based and model-free versions of the method

To illustrate the versatility of the idea, we will show an implementation of two versions of it in a simple case: the first that we will call model-based, refers to the case in which the functional form of the dynamic system and its symbolic dynamics are known, the second , which we will refer to as model-free, represents the case in which only data, resulting from the observation of the system, is available.

3.1 Model-based control

To illustrate our strategy, let us start considering the Logistic map, $f(s_{t})=4s_{t}(1-s_{t})$. It has an unstable fixed point at $s^{*}=3/4$ that we will consider as the target. The inverse function of $f$ can be written using the symbolic dynamic:

[TABLE]

The controlled system can be written as $s_{t+1}=4s_{t}(1-s_{t})+u_{t}$, with $u_{t}=u(x_{t})$. In this case, since the symbolic dynamic is known, the model-based and data-driven control schemes are showed. Figures 1 and 2, shows the results of the control scheme application using $500$ data points with $l=6$, $k=4$ and $q=16$ in the case of model-based control and $l=6$ and $q=16$ in the case of data-driven control.

It is worth to note that in general, if all possible sequences of symbolic states are considered as training data generation, the quality of the performance of the strategy could be affected, because this number of pairs does not necessarily coincide with the admissible ones. Now, if the dynamics are available, it is possible to determine the admissible symbolic sequences using the strategy proposed in [4]. The study of this aspect of the problem is not the subject of this article, but will be addressed in future works.

In this work, the strategy of control depend on the methodology to approximate $u$ but we believe that any strategy, to interpolate the response of the system to states that are not present in the their training set, can solve the control problem.

3.2 Model-free control

Since it is not always possible to have the symbolic dynamics of the system, in the following we shows how it is possible to implement our idea using data from a typical orbit of the system. Thus a data-driven version of the strategy is easily designed if we replace the symbolic dynamics knowledge with an algorithm that identifies, in the time series, a set of sequences of states converging to $s^{*}$.

At this point, we could use a numerical method of detection of homoclinic orbits as in references [3, 5, 17] and find the points which we will include in the target. However, we will opt for a simpler strategy: we start identifying, in the time series, the set of the nearest states to $s^{*}$, $\{s_{m_{1}},s_{m_{2}},\dots,s_{m_{q}}\}$, so that $s^{i}_{0}=s_{m_{i}-l}$, with $i=1,\cdots,q$ and $l$ is a integer identifying the $l$-th predecessor of $s_{m_{i}}$. Thus the set $S_{0}=\{s_{0}^{i}\}_{i=1}^{q}$ is constructed and given that set the control strategy is the same that in previous section.

In the following section the performance of the methodology is shown. As in the examples that follow it is not so simple to obtain the symbolic representation, we will only show results of the model-free approach.

The scheme can be summarized as the following list of tasks:

Give a time series $\{s_{i}\}$, chose: $\sigma^{*}$, $q$, $k$ and $l$ 2. 2.

Identify $q$ states $S_{0}=\{s^{1}_{0},s^{2}_{0},\cdots,s^{q}_{0}\}$, such that, with $s^{*}+\delta=f^{k}_{r}(s_{0})$.

• •

If the symbolic dynamic is known:

• i.

Identify a set of admissible states $\Sigma_{i}=\sigma^{i}_{1},\sigma^{i}_{2},\cdots,\sigma^{i}_{k},\underbrace{\sigma^{*},\cdots,\sigma^{*}}_{(l-k)~{}times}$,

• ii.

$s^{i}_{0}=f_{r}^{-l}(s_{R_{i}},\Sigma_{i})$ with $s_{R_{i}}$ a random real number.

• •

Else:

• i.

Identify a set of the nearest states to $s^{*}$, $\{s_{m_{1}},s_{m_{2}},\dots,s_{m_{q}}\}$,

• ii.

$s^{i}_{0}=s_{m_{i}-l}$. 3. 3.

Construct the training set $\{(s_{t},u_{t})_{i}\}_{i=1}^{N}$, with $u_{t}=P(S_{0},s_{t})-s_{t}$. 4. 4.

Train some Machine Learning Regression scheme, to obtain, $u_{t}=u(s_{t})$.

4 Numerical Results

In order to show the performance of the proposed scheme, in a more general context, we use as examples: the Henon map, the Lorenz system and a network of coupled discrete maps. These represent two, three and $n$-dimensional systems with discrete time and continuous time. Each of these systems has a homoclinic orbit, as can be seen in [28, 18]

4.1 Henon map

The Henon map is used as canonical example of bi-dimensional chaotic system and given by a difference equation system. It has an unstable fixed point at $(x^{*},y^{*})=(0.6313,0.1894)$ and the controlled system is given by:

[TABLE]

Here $u_{t}=u(x_{t},y_{t})$ and $10^{3}$ data points are used to to train the Gaussian Process with parameters, $l=9$ and $q=64$.

4.2 Lorenz system

The Lorenz system is the typical example of the continuous chaotic dynamical system represented by the set of ordinary differential equations. With parameters $\sigma=10$, $\rho=28$ and $\beta=8/3$, it has unstable fixed points at $(\pm\sqrt{\beta(\rho-1)},\pm\sqrt{\beta(\rho-1)},\rho-1)$.

The controlled version of the system, in this case is:

[TABLE]

Here $u(t)=u(x(t),y(t),z(t))$ and $1.2\times 10^{4}$ data points are used to to train the Random Forest Regression, with parameters, $l=10$ and $q=128$.

4.3 A coupled map network

The systems composed of parts that have some degree of autonomy and interact, are frequent in innumerable situations and in areas from basic to applied science, so they must be an obligatory example of control object. Here, without loss of generality we use the logistic map, $f(s^{i}_{t})=rs^{i}_{t}(1-s^{i}_{t})$, as the nonlinear local dynamic to construct a generic complex dynamical network consisting of $M$ identical nodes,

[TABLE]

where $L=(L_{ij})_{ij=1}^{M}$ is the usual Laplacian matrix with diagonal entries $L_{ii}=k_{i}$ and $k_{i}$ the out degree of the node $i$. Thus, the network (11) can be represented as

[TABLE]

where ${\bf s}_{t}$ are vectors, ${\bf f}({\bf s}_{t})=(f(s^{1}_{t}),f(s^{2}_{t}),\cdots,f(s^{M}_{t}))$ a vector valued function, ${\bf C}=(L_{ij}/k_{i}))_{ij=1}^{N}$ is a matrix with real eigenvalues[32] and $\epsilon$ a real parameter. Thus, the network can be represented as ${\bf s}_{t+1}={\bf A}{\bf f}({\bf s}_{t})$, with ${\bf A}={\bf I}-\epsilon~{}{\bf C}$.

This system have a unstable and homogeneous fixed point given by ${\bf s}^{*}=(s_{1}^{*},s_{2}^{*},\dots,s_{N}^{*})$. In this example $M=5$, the coupling parameter $\epsilon=0.1$ and

[TABLE]

Although the symbolic dynamics of this network can be constructed in a way similar to the one-dimensional case [22], here we will only present results for the data-dependent control strategy and although it is outside the our objectives the optimization of control sites we believe that this idea can be adapted to other contexts like the presented in [9].

The following figure shows evolution of the controlled system a a function of the time.

In this example, the controlled system is written as:

[TABLE]

where $u^{i}_{t}=u^{i}({\bf s}_{t})$.

Figures 6 and 7 shows the structure of the network used as example in this article and the evolution of the controlled network, respectively. The control functions $u^{i}_{t}$ was generated by training the machine learning regression with $2.5\times 10^{5}$ data points with parameters $l=7$ and $q=8\times 10^{3}$.

In spite of the results shown in Figure 7, refer to a network with the topology of the network in Figure 6, the performance of the control scheme is similar for many networks with the same number of units. Numerical experiments suggest that, as this number increases the amount of data needed to obtain good approximations of $u$ quickly grows so it would be necessary, in the case of large networks, to combine the control strategy with a numerical methodology to approximate homoclinic orbits. This is out of the scope of this work and it is currently under investigation.

5 Final remarks

A methodology for control of complex systems based on the estimation of the control function, from the observation of the system and using a machine learning technique, was presented. Although in this particular case we have used the regression method based on Gaussian processes or Random Forest, the methodology allows us to use the regression technique that best suits the particular problem.

Despite the existence of other control methods that use artificial intelligence strategi es, see for example [19] and references therein, th e proposed method is, as far as we know, novel in the sense that it directly approximate a nonlinear control signal from an rough estimate of homoclinic points and does so with very little computational effort. In addition, results suggest that the Proposed technique could be useful in the design of control or synchronization methodologies of complex systems like as [1].

Several aspects, although outside the scope of this work, remain as pending tasks. Among these aspects stand out: the optimization of the parameters ($q$, $k$ and $l$), selection of the sites to disturb or the design of a general criterion for the selection of regression strategy according to the number of data and dimension of the system.

Finally, the aim of the work is to draw attention to two ideas: i) it is possible to convert the control target into a set of targets, which in some cases would simplify the control task, ii) it is possible to approximate the control $u(s_{t})$ (not the system) using machine learning techniques. The results suggest that it is possible to implement both ideas into a simple control algorithm.