**1. Introduction**

Experts in the field of control are aware that the solution of the optimal control problem in the classical statement [1] cannot be implemented directly in a real control object, even in the presence of a sufficiently accurate mathematical model of the control object. The main reason for it is that the solution to the classical optimal control problem is a time function independent from coordinates of the state space vector. Therefore, the use of the obtained solution in the control object results in an open-loop control system that is sensitive even to small disturbances. The solution of the optimal control problem is only needed to obtain the optimal program trajectory in the state space. To provide a movement of a control object along this optimal trajectory, it is necessary to develop an additional control system with feedback control. It should be noted here that the development of an additional control system changes the mathematical model of the object that was used in the optimal control problem. Thus, the optimal trajectory, obtained for a control object, when solving an optimal control problem, is no longer optimal for an object with a control system that provides movement along an optimal trajectory.

There are several ways to solve the optimal control problem and provide the movement of an object along the optimal trajectory in the state space. In practice, engineers first make these robots stable relative to some point in space, and then control the movement of the robots by repositioning these stable points [2,3]. Usually, they set these points along a given trajectory to allow robots to move along that trajectory. They do not solve the optimal problem, because the control object stops near a stable point, and the task is to follow these

**Citation:** Diveev, A.; Sofronova, E.; Konyrbaev, N.; Bexeitova, A. Stabilization of Movement along an Optimal Trajectory and Its Solution. *Eng. Proc.* **2023**, *33*, 12. https:// doi.org/10.3390/engproc2023033012

Academic Editors: Ivan Zelinka, Arutun Avetisyan and Alexander Ilin

Published: 8 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

points subsequently. To ensure stability in the state space, engineers usually apply linear P-, PI- or PID-regulators in a feedback loop.

To provide practical feasibility of the solution, it is possible to solve the general control synthesis problem instead of the optimal control problem. For this purpose, it is necessary to change one initial state into a domain of initial states and to look for a control function of the state space vector. The search for a structure of function of many arguments is a difficult task. Here, it is possible to use methods of symbolic regression [4,5]. This perspective direction allows the search of mathematical expressions in coded form by means of a special evolutionary genetic algorithm. As a result, we obtain rather complex mathematical expressions, but they can be directly applied in a real control object.

Another way to solve the optimal control problem that can be implemented in a real object is to use a synthesized control specially developed for this purpose [6]. According to this approach, firstly, the control synthesis problem is solved, and a control object becomes stable in the state space relative to some equilibrium point. After that, the positions of this stable equilibrium point are defined by the quality criterion of the optimal control problem so that, if these positions are switched over equal time intervals, a control object every time trends to the stable equilibrium point and moves along the optimal trajectory. As a result, it achieves the given terminal state with the optimal value of the given quality criterion.

In this work, the third way of solution is considered, when after solving the optimal control problem a stabilization system for motion along the optimal trajectory is obtained. This problem is a control system synthesis for tracking the given trajectory. According to this approach, a control object is included in the control system as a reference model for the trajectory generation needed for tracking. Unlike some known works on tracking [7–10], here the optimal trajectories obtained by solving the optimal control problem are considered and a tracking system synthesis problem is solved by means of symbolic regression. In this work, mathematical statements of the extended optimal control problem [11] and the tracking system synthesis problem are presented. In the experimental part, the optimal control problem for a group of four quadcopters is presented.

#### **2. The Extended Optimal Control Problem**

Consider the optimal control problem with additional requirements to the optimal trajectory.

The mathematical model of control object is given

$$
\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{u}),
\tag{1}
$$

where **<sup>x</sup>** is a state-space vector of control object, **<sup>x</sup>** <sup>∈</sup> <sup>R</sup>*n*, **<sup>u</sup>** is a control vector, **<sup>u</sup>** <sup>∈</sup> <sup>U</sup> <sup>⊆</sup> <sup>R</sup>*m*, U is a compact set.

A compact set U determines restrictions on the control and is often presented in the form of algebraic inequalities

$$u\_i^- \le u\_i \le u\_i^+, \ i = 1, \dots, m,\tag{2}$$

where *ui* is a component of the control vector, **u** = [*u*<sup>1</sup> ... *um*] *<sup>T</sup>*, *u*<sup>−</sup> *<sup>i</sup>* , *<sup>u</sup>*<sup>+</sup> *<sup>i</sup>* are the given values, *i* = 1, . . . , *m*.

The initial state is given:

$$\mathbf{x}(0) = \mathbf{x}^0 \in \mathbb{R}^n. \tag{3}$$

The terminal state is given:

$$\mathbf{x}(t\_f) = \mathbf{x}^f \in \mathbb{R}^n,\tag{4}$$

where *tf* is a terminal time not given but limited, *tf* ≤ *t* <sup>+</sup>, *tf* is defined on achievement of the terminal state, *t* <sup>+</sup> is the given limit time of control process. If the control object does not reach the terminal state in *t* <sup>+</sup>, then it is considered that the control object will never reach the terminal state (4).

The quality criterion is given

$$J = \int\_0^{t\_f} f\_0(\mathbf{x}, \mathbf{u})dt \to \min\_{\mathbf{u} \in \mathcal{U}}.\tag{5}$$

It is necessary to find a control function in the following form

$$\mathbf{u} = \mathbf{v}(t), \ t \in (0; t\_f). \tag{6}$$

If we replace a control vector **u** by the found control function **v**(*t*) in the right part of ODE system (1), then the obtained ODE system

$$
\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, \mathbf{v}(t)) \tag{7}
$$

will have a particular solution that reaches the given terminal state (4) from the given initial state (3) with the optimal value of the quality criterion (5).

Let **v**∗(*t*) be the optimal control function. In the second stage, consider the following ODE system

$$\begin{array}{rcl} \dot{\mathbf{x}}^{\*} &=& \mathbf{f}(\mathbf{x}^{\*}, \mathbf{v}^{\*}(t)),\\ \dot{\mathbf{x}}^{\*} &=& \mathbf{f}(\mathbf{x}, \mathbf{u}). \end{array} \tag{8}$$

In (8), the first subsystem is the reference model that generates optimal trajectory. For the reference model, the initial state (3) is given.

For the second subsystem, a domain of initial states is given

$$
\chi\_0 \subseteq \mathbb{R}^n,\tag{9}
$$

where X0 is a compact set.

It is necessary to find a control function in the following form

$$\mathbf{u} = \mathbf{h}(\mathbf{x}^\* - \mathbf{x}) \in \mathbb{U},\tag{10}$$

where x∗ is the given optimal trajectory or the particular solution of the reference system in (8) from the given initial state (3).

The control function should minimize the following quality criterion

$$J\_1 = \int \cdots \int \max\_t \| \mathbf{x}^\*(t) - \mathbf{x}(t, \mathbf{y}) \| dt \to \min\_{\mathbf{h}(\mathbf{x}^\* - \mathbf{x}) \in \mathcal{U}} \tag{11}$$

where **y** ∈ X0.

#### **3. Symbolic Regression**

To solve the synthesis problem and find the control function (10), symbolic regression is used.

Symbolic regression is a perspective computational method that allows to find mathematical expressions fully automatically. If someday humanity manages to create artificial intelligence, that is, write a program that will reflect and form the goals of the calculation, then this will be performed not manually but automatically. Symbolic regression is such an approach. It allows finding the structure and parameters of mathematical expressions automatically in a coded form. Previously, when the problem of finding a mathematical expression arose, for example, for deriving the physical laws, humans wrote mathematical expressions with accuracy to values of some constant parameters. Only for finding these parameters on experimental data, numerical methods were used, for example, the least-squares method. Note, that any artificial neural network is a function with a given structure and a big number of unknown parameters. The search of these parameters by optimization algorithms is called machine learning.

In this work, a network operator method is used to solve the tracking system synthesis problem and to find the control function as a function of a state space vector.

To code mathematical expressions, symbolic regression uses an alphabet of elementary functions. The network operator uses an alphabet of unary and binary functions and codes a mathematical expression in the form of an oriented computational graph. Consider an example of coding a mathematical expression by the network operator.

Let us be given a mathematical expression

$$y = q\_1 \sin(\mathbf{x}\_1 + q\_2 \cos(\mathbf{x}\_1 - \mathbf{x}\_2)).\tag{12}$$

The expression has two parameters *q*1, *q*<sup>2</sup> and two variables *x*1, *x*2. The alphabet of the following elementary functions is enough to present this mathematical expression.

A set of unary functions

$$F\_1 = \{ f\_{1,1}(z) = z, f\_{1,2}(z) = -z, f\_{1,3}(z) = \sin(z), f\_{1,4}(z) = \cos(z) \},\tag{13}$$

A set of binary functions

$$F\_2 = \{f\_{2,1}(z\_1, z\_2) = z\_1 + z\_2, f\_{2,2}(z\_1, z\_2) = z\_1 z\_2\}.\tag{14}$$

Elementary function designations contain two indexes. The first index is the number of arguments of the function, and the second index is its number in the set. Using elementary functions, expression (12) can be written as

$$y = f\_{2,2}(q\_{1\prime}, f\_{1,3}(f\_{2,1}(\mathbf{x}\_1, f\_{2,2}(q\_2, f\_{1,4}(f\_{2,1}(\mathbf{x}\_1, f\_{1,2}(\mathbf{x}\_2)))))))).\tag{15}$$

In the network operator method, functions with two arguments *f*2,*i*(*z*1, *z*2) must be commutative, associative and have a unit element, *ei*. If one of the arguments is equal to a unit element, then the result of the calculations of the function is equal to the second argument, *f*2,*i*(*ei*, *z*2) = *z*2. For the addition function, the unit element is zero, *e*<sup>1</sup> = 0, and for multiplication, the unit element is one, *e*<sup>2</sup> = 1.

Figure 1 presents the network operator for the mathematical expression (15). Arguments of the mathematical expression are placed in the source nodes of the graph, the numbers of the binary functions are placed in other nodes of the graph and the numbers of the unary functions are located over arcs. To present the network operator in the PC memory, all nodes must be numerated so that the node number from which the arc outputs is less than the node number where the arc inputs. It can be done if the graph has no loops. Node numbers are pointed in the upper parts of the nodes.

**Figure 1.** The network operator of the mathematical expression.

The network operator is presented in the PC memory in the form of an upper triangular integer network operator matrix. The network operator matrix for the considered mathematical expression is

$$\mathbf{Y} = \begin{bmatrix} 0 & 0 & 0 & 0 & 1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 4 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 3 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 2 \end{bmatrix} \tag{16}$$

In the network operator matrix, the zero elements on the main diagonal point to the lines corresponding to the source nodes. The non-zero elements of the main diagonal are the binary function numbers. Non-zero elements above the main diagonal are the unary function numbers.

In order to calculate the mathematical expression by the network operator matrix, it is necessary to set a nodes vector of a dimension equal to the number of nodes of the network operator graph

$$\mathbf{z} = [z\_1 \dots z\_L]^T,\tag{17}$$

where *L* × *L* is a dimension of the network operator matrix.

Firstly, it is necessary to initialize the nodes vector. The components of the nodes vector corresponding to the source nodes are equal to the arguments of the mathematical expression, and the remaining components are equal to the unit elements of the corresponding binary functions.

Further, the values of the components of the nodes vector change after passing each line of the matrix and finding a non-zero non-diagonal element in it, as

$$z\_{i\_j}^{(i)} \leftarrow \begin{cases} f\_{2, \#\_{\left[i,j\right]}}(z\_j^{(i-1)}, f\_{1, \#\_{\left[i,j\right]}})(z\_i^{(i-1)}), \text{if } \psi\_{i,j} \neq 0\\ z\_j^{(i-1)} \end{cases}, i = 1, \dots, L-1, \; j = i+1, \dots, L,\tag{18}$$

where *ψi*,*<sup>j</sup>* is an element of the network operator matrix **Ψ** = [*ψi*,*j*], *i*, *j* = 1, . . . , *L*.

For the network operator matrix (16), a nodes vector has the following changes

$$\begin{array}{rcl} \mathbf{z}^{(0)} &=& [\mathbf{x}\_{1} \ x\_{2} q\_{1} q\_{2} \ 0 \ 1 \ 0 \ 1]^{T} \\ \mathbf{z}^{(1)} &=& [\mathbf{x}\_{1} \ x\_{2} q\_{1} q\_{2} \ \mathbf{x}\_{1} \ 1 \ \mathbf{x}\_{1} \ 1]^{T} \\ \mathbf{z}^{(2)} &=& [\mathbf{x}\_{1} \ x\_{2} q\_{1} q\_{2} \ \mathbf{x}\_{1} - \mathbf{x}\_{2} \ 1 \ \mathbf{x}\_{1} \ 1]^{T} \\ \mathbf{z}^{(3)} &=& [\mathbf{x}\_{1} \ x\_{2} q\_{1} q\_{2} \ \mathbf{x}\_{1} - \mathbf{x}\_{2} \ 1 \ \mathbf{x}\_{1} \ q\_{1}]^{T} \\ \mathbf{z}^{(4)} &=& [\mathbf{x}\_{1} \ x\_{2} q\_{1} q\_{2} \ \mathbf{x}\_{1} - \mathbf{x}\_{2} q\_{2} \ \mathbf{x}\_{1} q\_{1}]^{T} \\ \mathbf{z}^{(5)} &=& [z\_{1}^{(4)} \ \dots z\_{5}^{(4)} \ q\_{2} \cos(\mathbf{x}\_{1} - \mathbf{x}\_{2}) \ \mathbf{x}\_{1} q\_{1}]^{T} \\ \mathbf{z}^{(6)} &=& [z\_{1}^{(5)} \ \dots z\_{6}^{(5)} \ \sin(\mathbf{x}\_{1} + q\_{2} \cos(\mathbf{x}\_{1} - \mathbf{x}\_{2})) \ q\_{1}]^{T} \\ \mathbf{z}^{(7)} &=& [z\_{1}^{(6)} \ \dots z\_{7}^{(6)} \ \$$

For more details on the network operator and the genetic algorithm that searches for the mathematical expression using the network operator, see the monographs [4,5].

#### **4. Computational Experiment**

Consider the optimal control problem for spatial movement of four similar quadcopters.

The mathematical model of control object is

$$\begin{array}{rcl} \dot{x}\_1^j &=& \dot{x}\_4^j, \\ \dot{x}\_2^j &=& \dot{x}\_5^j, \\ \dot{x}\_3^j &=& \dot{x}\_6^j, \\ \dot{x}\_4^j &=& u\_4^j (\sin(u\_3^j)\cos(u\_2^j)\cos(u\_1^j) + \sin(u\_1^j)\sin(u\_2^j)), \\ \dot{x}\_5^j &=& u\_4^j \cos(u\_3^j)\cos(u\_1^j) - \xi\_{\nu'} \\ \dot{x}\_6^j &=& u\_4^j (\cos(u\_2^j)\sin(u\_1^j) - \cos(u\_1^j)\sin(u\_2^j)\sin(u\_3^j)), \end{array} \tag{20}$$

where *j* = 1, . . . , 4, *gc* = 9.80665.

The control is constrained

$$\begin{array}{rclrclrcl} -\pi/12 & \leqslant & u\_1^j & \leqslant & \pi/12, \\ -\pi & \leqslant & u\_2^j & \leqslant & \pi, \\ -\pi/12 & \leqslant & u\_3^j & \leqslant & \pi/12, \\ 0 & \leqslant & u\_4^j & \leqslant & 12, \quad j = 1, \ldots, 4. \end{array} \tag{21}$$

For the system (20), the initial states are given:

$$\begin{array}{rcl} \mathbf{x}^{1,0} &=& [0 \ 5 \ 0 \ 0 \ 0 \ 0]^T, & \mathbf{x}^{2,0} &=& [10 \ 5 \ 0 \ 0 \ 0 \ 0]^T, \\ \mathbf{x}^{3,0} &=& [0 \ 5 \ 10 \ 0 \ 0 \ 0]^T, & \mathbf{x}^{1,0} &=& [10 \ 5 \ 10 \ 0 \ 0 \ 0]^T. \end{array} \tag{22}$$

The terminal states are given

$$\begin{array}{rcl} \mathbf{x}^{1,f} & = & [10 \ 5 \ 10 \ 0 \ 0 \ 0]^T, & \mathbf{x}^{2,f} & = & [0 \ 5 \ 10 \ 0 \ 0 \ 0]^T, \\ \mathbf{x}^{3,f} & = & [10 \ 5 \ 0 \ 0 \ 0 \ 0]^T, & \mathbf{x}^{1,f} & = & [0 \ 5 \ 0 \ 0 \ 0 \ 0]^T. \end{array} \tag{23}$$

The static phase constraints are given

$$\varphi\_i(\mathbf{x}^j) = r\_i - \sqrt{(\mathbf{x}\_i - \mathbf{x}\_1^j)^2 + (z\_i - \mathbf{x}\_3^j)^2} \tag{24}$$

where *i* = 1, ... , 4, *j* = 1, ... , 4, *r*<sup>1</sup> = *r*<sup>2</sup> = *r*<sup>3</sup> = *r*<sup>4</sup> = 1.5, *x*<sup>1</sup> = 1.5, *z*<sup>1</sup> = 2.5, *x*<sup>2</sup> = 1.5, *z*<sup>2</sup> = 7.5, *x*<sup>3</sup> = 8.5, *z*<sup>3</sup> = 2.5, *x*<sup>4</sup> = 8.5, *z*<sup>4</sup> = 7.5.

The dynamic phase constraints are given

$$\chi(\mathbf{x}^{j\_1}, \mathbf{x}^{j\_2}) = s - \sqrt{\sum\_{i=1}^{3} (\mathbf{x}\_i^{j\_1} - \mathbf{x}\_i^{j\_2})^2} \ll 0,\tag{25}$$

where *j*1, *j*<sup>2</sup> ∈ {1, 2, 3, 4}, *j*<sup>1</sup> = *j*2, *s* = 1.5.

Initially, the optimal control problem is solved by the direct approach. For this purpose, the time axis is divided into equal time intervals. In each interval, a control function is found as a piecewise-linear time function. It is necessary to find the values of the constant parameters on the boundaries of the time intervals. Taking into account the constraints on control, the control function is

$$u\_i^j = \begin{cases} u\_i^{j,+}, \text{if } \vec{u}\_i^j > u\_i^{j,+} \\ u\_i^{j,-}, \text{if } \vec{u}\_i^j < u\_i^{j,-} \\ \vec{u}\_{i'}^j, \text{otherwise} \end{cases},\tag{26}$$

where

$$\pi\_i^j = q\_{i + (k-1)m} + (q\_{i + km} - q\_{i + (k-1)m}) \frac{t - (k-1)\Delta t}{\Delta t}, \quad (k-1)\Delta t \lessapprox t \lessapprox k\Delta t,\tag{27}$$

*i* = 1, ... , *m*, *m* is a dimension of the control vector, *m* = 4, *k* = 1, ... , *N*, *N* is the number of intervals.

It is necessary to find the control functions for each control object as (26). The time limit of the control process is *t* <sup>+</sup> = 5.6. The time interval is Δ*t* = 0.4. To solve the optimal control problem, it is necessary to find 240 parameters, 4 · *m* ·(*N* + 1) = 4 · 4 · 5.6 0.4 + 1 = 4 · 4 · 15 = 240.

To solve the optimal control problem by direct approach, the evolutionary hybrid algorithm is used [6,12]. The algorithm includes evolutionary transformations of three evolutionary algorithms, the genetic algorithm (GA) [13], particle swarm optimization (PSO) algorithm [14] and grey wolf optimizer (GWO) algorithm [15].

Figure 2 shows projections of the found optimal trajectories on the horizontal plane. Here, the solid line is for the trajectory of the first robot, dashed line—the second robot, dotdashed line—the third robot, dots—the forth robot, and circles are the static phase constraints.

**Figure 2.** Optimal trajectories on the horizontal plane.

At the second stage, the tracking system synthesis problem is solved. To solve the synthesis problem, we used the network operator method.

The domain of initial states (9) was replaced by the set of points

$$\mathbf{x}\_{i}(0) = \mathbf{x}\_{i}^{j,0} \pm 0.2, \quad i = 1, 2, 3, \quad j = 1, \ldots, 4. \tag{28}$$

In the quality criterion (11), an integral over domain was replaced by the sum of all 33 = 27 initial states.

The following solution was found by the network operator method

$$u\_i^j = \begin{cases} u\_i^{j,+} \text{, if } \hat{u}\_i^j > u\_i^{j,+} \\ u\_i^{j,-} \text{, if } \hat{u}\_i^j < u\_i^{j,-} \\ u\_{i'}^j \text{ otherwise} \end{cases} \tag{29}$$

where

$$\mathcal{U}\_1^{\bar{j}} = \mu(A),\tag{30}$$

$$
\hat{\mu}\_2^j = \mu(A) - \mu^3(A), \tag{31}
$$

$$\mathfrak{u}\_3^j = \rho\_{17}(\mathfrak{u}\_2^j) + \rho\_{19}(B + \mu(A)) + \rho\_{17}(\mathbb{C}),\tag{32}$$

*u*ˆ *j* <sup>4</sup> = *u*ˆ *j* <sup>3</sup> + ln(*u*ˆ *j* <sup>2</sup>) + sgn(*<sup>B</sup>* <sup>+</sup> *<sup>μ</sup>*(*A*))|*<sup>B</sup>* <sup>+</sup> *<sup>μ</sup>*(*A*)<sup>|</sup> <sup>+</sup> *<sup>ρ</sup>*19(*B*)+ arctan(*D*) + sgn(*E*) + arctan(*F*) + exp(*q*2(*xj*,<sup>∗</sup> <sup>2</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>2</sup>)) + <sup>√</sup>*q*1, (33) *<sup>A</sup>* <sup>=</sup> *<sup>q</sup>*3(*xj*,<sup>∗</sup> <sup>3</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>3</sup>) + *<sup>q</sup>*6(*xj*,<sup>∗</sup> <sup>6</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* 6), *B* = *G* + tanh(*E* + sin(*C*) <sup>√</sup><sup>3</sup> *<sup>F</sup>*) + exp(*H*) <sup>−</sup> *<sup>C</sup>*, *<sup>C</sup>* <sup>=</sup> *<sup>q</sup>*1(*xj*,<sup>∗</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>1</sup>) + *<sup>q</sup>*4(*xj*,<sup>∗</sup> <sup>4</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* 4) *D* = *E* + sin(*C*) <sup>√</sup><sup>3</sup> *<sup>F</sup>* <sup>+</sup> tanh(*E*) + *<sup>ρ</sup>*18(*V*), *E* = *F* + arctan(*A*) + arctan(*H*) − *V*, *F* = *H* + sgn(*xj*,<sup>∗</sup> <sup>5</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* 5)+(*xj*,<sup>∗</sup> <sup>2</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>2</sup>)3, *<sup>G</sup>* <sup>=</sup> exp(*D*) + cos(*q*6(*xj*,<sup>∗</sup> <sup>6</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>6</sup>)) + sgn(*E* + sin(*C*) <sup>√</sup><sup>3</sup> *<sup>F</sup>*) |*E* + sin(*C*) <sup>√</sup><sup>3</sup> *<sup>F</sup>*<sup>|</sup> <sup>+</sup> *<sup>E</sup>*−<sup>1</sup> <sup>−</sup> *<sup>q</sup>*6, *<sup>H</sup>* <sup>=</sup> *<sup>ρ</sup>*17(*A*) + *<sup>V</sup>*<sup>3</sup> <sup>+</sup> *<sup>C</sup>* <sup>+</sup> *<sup>ϑ</sup>*(*q*5(*xj*,<sup>∗</sup> <sup>5</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>5</sup>)) + *<sup>x</sup>j*,<sup>∗</sup> <sup>5</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* 5, *<sup>V</sup>* <sup>=</sup> sin(*q*6(*xj*,<sup>∗</sup> <sup>6</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>6</sup>)) + *<sup>q</sup>*5(*xj*,<sup>∗</sup> <sup>5</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* <sup>5</sup>) + cos(*q*1) + *<sup>ϑ</sup>*(*xj*,<sup>∗</sup> <sup>2</sup> <sup>−</sup> *<sup>x</sup><sup>j</sup>* 2),

$$\theta(\mathfrak{a}) = \begin{cases} 1, \text{if } \mathfrak{a} > 0 \\ 0, \text{otherwise} \end{cases}, \quad \mu(\mathfrak{a}) = \begin{cases} \mathfrak{a}, \text{if } |\mathfrak{a}| < 1 \\ \text{sgn}(\mathfrak{a}), \text{otherwise} \end{cases}$$

*ρ*17(*α*) = sgn(*α*)ln(|*α*| +1), *ρ*18(*α*) = sgn(*α*)(exp(|*α*|) −1), *ρ*19(*α*) = sgn(*α*) exp(−|*α*|), *q*<sup>1</sup> = 11.52295, *q*<sup>2</sup> = 7.54395, *q*<sup>3</sup> = 15.45825, *q*<sup>4</sup> = 12.39771, *q*<sup>5</sup> = 15.94922, *q*<sup>6</sup> = 10.84229, **x***j*,<sup>∗</sup> = [*xj*,<sup>∗</sup> <sup>1</sup> ... *<sup>x</sup>j*,<sup>∗</sup> 6 ] *<sup>T</sup>* is a program trajectory, generated by the reference model.

Figure 3 shows the trajectories from eight random initial states of given domain (28). Here, blue lines mark trajectories from (3) and black lines mark trajectories from (28). As it can be seen from the experiment, the obtained control system (29) for optimal trajectory tracking guarantees the presence of an attracting neighborhood for the optimal trajectory in the extended optimal control problem.

**Figure 3.** The optimal trajectory and disturbances trajectories from eight initial states.

#### **5. Results**

The paper presents the statements of the optimal control problem and the stabilization system synthesis with respect to the optimal program trajectory. To solve the synthesis problem, machine learning by the network operator method was used. A brief description of the network operator method is presented. A computational example is given for a group of four quadcopters. To analyze the quality of the obtained solutions, we simulated the synthesized stabilization system under perturbations of the initial state. The experiments showed the high quality of the stabilization system.

#### **6. Discussion**

The implementation of the control system is related to the solution of the control synthesis problem. As a result of solving the synthesis problem, we obtain a feasible feedback control. It is a complex problem that does not have universal computational methods. In this paper, we applied an approach based on the machine learning by the symbolic regression method. The method is universal, but requires lots of calculations. It is necessary to continue the study of this approach for other complex control problems.

**Author Contributions:** Conceptualization, A.D.; methodology, A.D.; software, A.D. and E.S.; validation, A.D., E.S., N.K. and A.B.; investigation, A.D. and E.S.; writing—original draft preparation, A.D., E.S. and N.K.; writing—review and editing, E.S.; visualization, A.B.; supervision, N.K.; project administration, N.K.; funding acquisition, N.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Science Committee of the Ministry of Education and Science of the Republic of Kazakhstan (Grant No. AP14869851).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Experimental data is available at: https://cloud.mail.ru/public/L51J/ oDPtP3Xhc (accessed on 7 June 2023).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
