Solving the Control Synthesis Problem Through Supervised Machine Learning of Symbolic Regression

Abstract

This paper considers the control synthesis problem and its solution using symbolic regression. Symbolic regression methods, which were previously called genetic programming methods, allow one to use a computer to find not only the parameters of a given regression function but also its structure. Unlike other works on solving the control synthesis problem using symbolic regression, the novelty of this paper is that for the first time this work employs a training dataset to address the problem of general control synthesis. Initially, the optimal control problem is solved from each point in a given set of initial states, resulting in a collection of control functions expressed as functions of time. A reference model is then integrated into the control object model, which generates optimal motion trajectories using the derived optimal control functions. The control synthesis problem is framed as an approximation task for all optimal trajectories, where the control function is sought as a function of the deviation of the object from the specified terminal state. The optimization criterion for solving the synthesis problem is the accuracy of the object’s movement along the optimal trajectory. The paper includes an example of solving the control synthesis problem for a mobile robot using a supervised machine learning method. A relatively new method of symbolic regression, the method of variational complete binary genetic programming, is studied and proposed for the solution of the control synthesis problem.

1. Introduction

The problem of control system synthesis is an important and quite challenging issue when designing automatic control systems for autonomous objects. The general problem of control synthesis was first stated in the 1960s by V.G. Boltyansky [1] as a problem of an optimal control search from any initial condition. At that time, the solution was searched analytically using Pontryagin’s Maximum Principle [2] as a general solution for differential equations of control objects and conjugate variables with optimal control in the right parts. Not all differential equations have general solutions.

The complexity of the control synthesis problem lies in the fact that it is necessary to find a multidimensional control function, the argument of which is the state vector of the control object. If this function is substituted into the right parts of the differential equations of the control object model, then the resulting system of equations will have the following property: any particular solution from the initial state of the state space domain will reach the terminal state with the optimal value of the quality criterion.

For linear autonomous systems, the synthesis of the control system involves designing negative feedback based on the state error [3]. Today, the control synthesis problem is typically solved manually by engineers and control specialists to implement the developed control system. Engineers analyze the problem and the control object and establish control channels, i.e., determine the necessary control actions to compensate for specific state errors. Then, to improve control quality, regulators are inserted into these channels [4,5,6]. Most commonly, linear proportional (P), proportional–integral (PI) and proportional–integral–derivative (PID) controllers are used, each with their own gain coefficients. Traditional PID control algorithms feature simple structures, ease of implementation, and, at the same time, fixed parameters and lack of adaptability to complex environments [7]. To eliminate drawbacks, the PID control is enhanced by other control technologies such as fuzzy logic [8,9], etc.

This work employs a numerical approach to the control synthesis problem based on machine learning using symbolic regression [10,11,12,13,14]. In mathematics, a regression method was previously used to approximate experimental data. The researcher studied experimental data and wrote down the mathematical expression of the function, which approximates these data up to some parameters. The parameter values were then found by the least squares method. This type of regression is subjective and only possible for functions with few arguments. Symbolic regression is a modern computational method that allows you to find not only parameters, but the structure of a function. For this purpose, an alphabet of elementary functions is defined and the rules for encoding mathematical expressions using this alphabet are stipulated. Further, a special genetic algorithm searches for the optimal expression according to a given criterion of a mathematical expression in the code space. In the genetic algorithm, the main operation of the crossover ensures that codes of new mathematical expressions appear. All symbolic regression methods differ in the form of coding.

It is important to note that the control function is often a multidimensional nonlinear function of the vector argument. In practice, artificial neural networks are rarely used to approximate this function, as it appears on the right-hand side of the model equations and significantly alters the dynamic properties of the object. Moreover, since a mathematical expression for this function is not known in advance, it is impossible to obtain a sufficient training dataset to train a neural network.

Since 2008, symbolic regression methods have been utilized to address the control synthesis problem [15]. The control synthesis problem includes finding a control function as a function of state. It is often solved by supplying a stability property of the control object relative to a point in the state space. Then, for non-complex mathematical models, analytical methods can be used, such as a backstepping integrator [16,17,18] or dynamic programming [19,20,21], that provide a numerical value of control for each state of the control object, but one cannot obtain a mathematical expression for optimal control function in the appropriate form.

Here, the control synthesis problem is solved by symbolic regression. These methods involve encoding mathematical expressions based on a selected alphabet of elementary functions and searching the space of codes for the optimal mathematical expression of the function according to a specified quality criterion using a specialized genetic algorithm.

Research on the application of symbolic regression methods for solving the control synthesis problem has shown that most symbolic regression methods cannot directly address this task. The crossover operation leads to significant changes in codes, ultimately generating entirely new potential solutions that do not retain the properties of the selected parent solutions. Instead of producing offspring solutions that inherit the characteristics of their parents, completely new solutions are generated; thus, the genetic crossover operation functions more like a generator of new possibilities, transforming the evolutionary genetic algorithm into a simple random search. For the successful application of symbolic regression in finding the structure of mathematical expressions, it is essential for the crossover operation to possess the inheritance property with specific characteristics. To ensure this property, a principle of small variations in the basic solution was developed [22], and it has been demonstrated that implementing this principle allows the crossover operation to exhibit inheritance properties. This enhances the likelihood of finding the optimal solution compared to a simple random search.

The application of the principle of small variations in the baseline solution to the problem of general control synthesis has led to the modification of established symbolic regression methods. Further research on solving the general control synthesis problem using symbolic regression has revealed that these methods produce solutions in the form of encoded mathematical expressions. However, these solutions can differ significantly not only in terms of the optimization criterion value but also in the quality of the solutions themselves.

Unlike previous studies [23], this work employs a training dataset to address the problem of general control synthesis. Initially, the optimal control problem is solved from each point in a given set of initial states, resulting in a collection of control functions expressed as functions of time. A reference model is then integrated into the control object model, which generates optimal motion trajectories using the derived optimal control functions. The control synthesis problem is framed as an approximation task for all optimal trajectories, where the control function is sought as a function of the deviation of the object from the specified terminal state. The optimization criterion for solving the synthesis problem is the accuracy of the object’s movement along the optimal trajectory. The work includes an example of solving the control synthesis problem for a mobile robot using a supervised machine learning method, namely original variational complete binary genetic programming.

2. The Problem of General Control Synthesis

In the general control synthesis problem, it is necessary to find a control function whose argument is a state vector. The replacement of the control vector on the right part of the model with a control function should change the properties of the system so that all its particular solutions from some state space domain must reach a terminal state with an optimal value of the quality criterion.

The mathematical model of the control object in the form of an ODE system is given as follows:

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

(1)

where $x$ is a control state vector, $x\in {ℝ}^n$, $u$ is a control vector, $u\in \mathrm{U}\subseteq {ℝ}^m$ and $\mathrm{U}$ is a compact set that often defines the constraints on the control.

$$u^{-}\le u\le u^{+},$$

(2)

where $u^{-}$ and $u^{+}$ are the lower and upper boundaries of the control vector, respectively.

The terminal state is given as

$$x(t_f)=x^f={[x_1^f\dots x_n^f]}^{\mathrm{T}},$$

(3)

where $t_f$ is a time of reaching the terminal state, $t_f\le t^{+}$ and $t^{+}$ is a given positive value.

The set of initial states is given as

$${\mathrm{X}}_0=\left\{x^{0,1},x^{0,2},\dots ,x^{0,k}\right\},$$

(4)

$$x^{0,j}={[x_1^{0,j}\dots x_n^{0,j}]}^{\mathrm{T}}.$$

(5)

$$J_0={\displaystyle \sum _{j=1}^K\left({\displaystyle \underset{0}{\overset{t_{f,j}}{\int }}{f}_0(x,u)dt+p_1{‖x^f-x(t_{f,j})‖}_2}\right)}\to \underset{u\in \mathrm{U}}{\min },$$

(6)

where $p_1$ is a penalty coefficient,

$$t_{f,j}=\left\{\begin{array}{l}t, \mathrm{if} t<t^{+} \mathrm{and} \left|\right|x^f-x(t)|{|}_2\le \epsilon \hfill \\ t^{+}, \mathrm{otherwise}\hfill \end{array}\right., j=1,\dots ,K,$$

(7)

where $\epsilon$ is a given small positive value.

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

$$u=h(x^f-x),$$

(8)

that satisfies the restrictions on the control (2) and supplies the minimum of the quality criterion (6).

Our approach includes solving the optimal control problem. Therefore, the statement of the optimal control problem is presented here.

Analytically, this problem, (1)–(8), is solved only for very simple models and criteria, since the solution is a multidimensional function of a multidimensional argument, and this function is almost always nonlinear. The universal numerical method became possible only after the development of symbolic regression methods. These methods allow the search for the structure and parameters of any function according to a given criterion. In this study, we use a symbolic regression method to approximate many solutions to optimal control problems.

The mathematical model of a control object in the form of an ODE system is given (1). The constraints on the control (2) are given. The initial state is given as

$$x(0)=x^0={[x_1^0\dots x_n^0]}^{\mathrm{T}}.$$

(9)

The terminal state (3) is given. The quality criterion is given as

$$J_1={\displaystyle \underset{0}{\overset{t_f}{\int }}{f}_0(x,u)dt}+p_1||x^f-x(t_f)||\to \underset{u\in \mathrm{U}}{\min },$$

(10)

where $t_f$ is a time of reaching the terminal state, and it is described by Equation (7) for one initial state, $K=1$.

It is necessary to find a control function as a function of time satisfying the constraints on the control (2).

$$u=v(t), 0\le t\le t_f\le t^{+},$$

(11)

If the found optimal control function (11) is placed in the right part of the mathematical model (1), then the ODE system,

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

(12)

will have the particular solution $x(t,x^0)$ from the initial state (9), which reaches the terminal state (3) with the optimal value of the given quality criterion (10).

Firstly, the optimal control problem is solved for each initial state from the given set (4) of initial states. As a result, a set of optimal control functions, parameterized by time, is obtained.

$${\mathrm{V}}^{*}=\{v^{*,1}(t),\dots ,v^{*,K}(t)\}.$$

(13)

Afterward, the reference model is integrated into the mathematical model of the control object.

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

(14)

$${\dot{x}}^{*}=f(x^{*},v^{*}(t)).$$

(15)

Any particular solution of the reference model (15) for an initial state from the given set of initial states (4), along with the corresponding optimal control function from the set of optimal control functions (13), results in the generation of the optimal trajectory in the state space.

$$x^{*}(t)=x^{*}(t,x^{0,j}), j=1,\dots ,K,$$

(16)

where $x^{0,j}\in {\mathrm{X}}_0$, $v^{*,j}(t)\in {\mathrm{V}}^{*}$ and $j\in \{1,\dots ,K\}$.

In the second stage, the control synthesis problem in (1)–(6) and (8) is addressed with the following quality criterion:

$$J_1={\displaystyle \sum _{j=1}^K\left({\displaystyle \underset{0}{\overset{t_{f,j}}{\int }}\left|\right|x^{*}(t,x^{0,j})-x(t,x^{0,j})|{|}_{2}dt}+\underset{t}{\max \left|\right|x^{*}(t,x^{0,j})-x(t,x^{0,j})|{|}_2}\right)}+$$

$$p_1{\displaystyle \sum _{j=1}^K\left(\left|\right|x^f-x(t_{f,j},x^{0,j}|{|}_2\right)\to \underset{h(x^f-x)\in \mathrm{U}}{\min }},$$

(17)

where $t_{f,j}$, $j\in \{1,\dots ,K\}$ is defined by Equation (7).

3. Symbolic Regression for Solving the Control Synthesis Problem

Symbolic regression encodes all mathematical expressions in the form of a special code and utilizes a genetic algorithm to search for the optimal mathematical expression. Let us introduce variational complete binary genetic programming (VCBGP).

Consider examples of coding a mathematical expression using this method. Consider the following mathematical expression:

$$y=\sqrt[3]{-{\displaystyle \frac{x_2}{q_1}}+\sqrt{{\displaystyle \frac{x_2^2}{q_1^2}}+{\displaystyle \frac{x_1^3}{q_2^3}}}}+\sqrt[3]{-{\displaystyle \frac{x_2}{q_1}}-\sqrt{{\displaystyle \frac{x_2^2}{q_1^2}}+{\displaystyle \frac{x_1^3}{q_2^3}}}}.$$

(18)

This is Cardano’s formula for solving the cubic equation

$$y^3+x_{1}y+x_2=0,$$

at $q_1=2$ and $q_2=3$.

To encode this mathematical expression, an alphabet of elementary functions is required. For the given mathematical expression (18), the following set of elementary functions is sufficient:

-

Functions with one argument:

$${\mathrm{F}}_1=\{f_{1,1}(z)=z, f_{1,2}(z)=-z, f_{1,3}(z)=\sqrt{z}, f_{1,4}(z)=z^2,$$

$$f_{1,5}(z)=z^{-1}, f_{1,6}(z)=\sqrt[3]{z}, f_{1,7}(z)=z^3\},$$

(19)

-

Functions with two arguments:

$${\mathrm{F}}_2=\{f_{2,1}(z_1,z_2)=z_1+z_2, f_{2,2}(z_1,z_2)=z_1{z}_2\}$$

(20)

-

Arguments of mathematical expression or functions without arguments:

$${\mathrm{F}}_0=\{f_{0,1}=x_1, f_{0,2}=x_2, f_{0,3}=q_1, f_{0,4}=q_2, f_{0,5}=q_3, f_{0,6}=1\}.$$

(21)

To construct a graph of VCBGP, it is first necessary to create a graph of binary genetic programming and then extend it to form the complete binary genetic programming. To augment the graph of binary genetic programming with the required number of levels, it is essential to insert an addition function with zero arguments or a multiplication function with unit arguments into the graph. Therefore, the set of arguments includes unit elements for both the addition and multiplication functions.

All elements of the alphabet have two subscripts: the first subscript denotes the number of arguments, while the second subscript indicates the element’s position in the set. The set of functions with one argument (19) includes the identity function $f_{1,1}(z)=z$. All functions with two arguments (20) are commutative, associative and have their own unit element.

Binary genetic programming encodes a mathematical expression as a computational tree. This method employs only functions with one and two arguments. The representation of a mathematical expression in the form of binary genetic programming is illustrated in Figure 1.

In Figure 1, the nodes of the tree display the numbers of the functions with two arguments, while the edges of the graph are labeled with the numbers of the functions with one argument. The arguments of the mathematical expression are represented at the last level of the graph or at the leaves of the tree.

Complete binary genetic programming differs from standard binary genetic programming in that it has the maximum possible number of nodes at all levels. At level $l$ in complete binary genetic programming, there should be $2^{l-1}$ nodes.

According to the graph of binary genetic programming, the number of levels required to encode the mathematical expression (18) should be at least $L=6$. In the binary genetic programming graph, at the third level, the first and third functions with two arguments do not contain a second argument. Figure 2 shows the completion of the binary genetic programming graph with the second argument for functions with two arguments.

Figure 2a presents a part of the graph from Figure 1, where a function with two arguments uses only one argument. Figure 2b shows the correction of the graph to represent it in the form of complete binary genetic programming.

In Figure 2 of the binary genetic programming graph, at the third level, the arguments of the functions at nodes 2 and 4 are the arguments of the mathematical expression itself, although these arguments should be located at the sixth level of the complete binary genetic programming graph.

Figure 3a presents a part of the binary genetic programming graph from Figure 1, which contains the arguments of the mathematical expression at the fourth level. Figure 3b shows the same part of the graph, but with two additional levels added on top.

The redundancy of nodes in the complete binary genetic programming graph, compared to standard binary genetic programming, is necessary to more efficiently represent the graph in computer memory. Since the number of elements at each level in the complete binary genetic programming graph is known, representing this graph in computer memory requires only an ordered set of integers, which indicates the function numbers at each level from left to right. By sequentially listing the function numbers from left to right, it becomes easy to determine the number of arguments for each function based on its number. In total, the code for complete binary genetic programming for a graph with $L$ levels must contain the following elements:

$$S={\displaystyle \sum _{i=1}^L{2}^i}=2+4+\cdots +2^L=2^{L+1}-2$$

(22)

Accordingly, the last $2^{L-1}$ elements represent the arguments of the mathematical expression.

Let us consider the code in complete binary genetic programming.

$$\mathrm{B}=(b_1,\dots ,b_M).$$

(23)

Element $b_i$, $i\in \{1,\dots ,M=2^{L+1}-2\}$, encodes the following elements of the mathematical expression:

• If $2^L<i\le 2^{L+1}$, then $b_i$ is the argument of the mathematical expression;
• If $2^k<i\le 2^{k+1}$ and $k\mathrm{mod}2=0$,$k<L$, then $b_i$ is the number of the function with two arguments;
• If $2^k<i\le 2^{k+1}$ and $k\mathrm{mod}2\ne 0$, $k<L$, then $b_i$ is the number of the function with one argument.

To determine the number of levels, we use the following relation:

$$L={\log }_2(M+2)-1,$$

(24)

The code for the complete binary genetic programming of the mathematical expression (18) has the following form:

$$\mathrm{B}=(1,1,/6,6,1,1,/1,2,2,2,1,1,1,1,/3,1,1,1,3,1,1,1,1,1,1,1,1,1,1,1,/$$

$$4,7,1,1,1,1,1,1,4,7,1,1,1,1,1,1,2,2,1,1,1,1,1,2,2,2,1,1,1,1,1,2,/$$

$$5,1,5,1,1,1,1,1,1,1,1,1,1,1,5,1,5,1,5,1,1,1,1,1,1,1,1,1,1,1,$$

$$3,2,4,1,5,5,5,5,5,5,5,5,5,5,3,2,3,2,4,1,5,5,5,5,5,5,5,5,5,5,3,2),$$

(25)

The layers in (25) are separated by slashes.

In search of the mathematical expression, the VCBGP employs a variational genetic algorithm. This algorithm utilizes the principle of small variations in the basic solution. According to this principle, the initial population encodes only one basic solution using the code of VCBGP. All other possible solutions are encoded by an ordered set of codes representing small variations.

In VCBGP, an integer vector of two components is used to encode a small variation,

$$w={[w_1 w_2]}^{\mathrm{T}},$$

(26)

where $w_1$ is the position number in code, and $w_2$ is a new value of the element in the pointed position.

In the population, each possible solution is represented as an ordered set of small variation vectors,

$${\mathrm{W}}_i=(w^{i,1},\dots ,w^{i,d}), i=1,\dots ,H,$$

(27)

where $d$ is the depth of variation; $H$ is the number of possible solutions in the population except the basic solution.

Any possible solution, represented by the VCBGP code ${\mathrm{B}}_i,$ is a result of applying a set of small variations to the basic solution,

$${\mathrm{B}}_i={\mathrm{W}}_i\circ {\mathrm{B}}_0=w^{i,d}\circ \dots \circ w^{i,1}\circ {\mathrm{B}}_0,$$

(28)

where ${\mathrm{B}}_0$ is a VCBGP code of the basic solution.

During the crossover operation, two possible solutions, represented as ordered sets of small variation vectors, are selected randomly.

$${\mathrm{W}}_{\alpha }=(w^{\alpha ,1},\dots ,w^{\alpha ,d}), {\mathrm{W}}_{\beta }=(w^{\beta ,1},\dots ,w^{\beta ,d}),$$

(29)

A crossover point is determined randomly, $c\in \{1,\dots ,d\}$.

Two new possible solutions are generated by exchanging the tails of the selected solutions after the crossover point.

$${\mathrm{W}}_{H+1}=(w^{\alpha ,1},\dots ,w^{\alpha ,c},w^{\beta ,c+1},\dots ,w^{\beta ,d}),$$

(30)

$${\mathrm{W}}_{H+2}=(w^{\beta ,1},\dots ,w^{\beta ,c},w^{\alpha ,c+1},\dots ,w^{\alpha ,d})$$

(31)

4. Computational Experiment

Consider the synthesis of a stabilization system for a wheeled mobile robot with differential drive. The mathematical model of the control object is composed of three ordinary differential equations,

$${\dot{x}}_1=0.5(u_1+u_2)\cos (x),$$

(32)

$${\dot{x}}_2=0.5(u_1+u_2)\sin (x_3),$$

(33)

$${\dot{x}}_3=0.5(u_1-u_2),$$

(34)

where $x={[x_1 x_2 x_3]}^{\mathrm{T}}$ is a state vector of the mobile robot, $u={[u_1 u_2]}^{\mathrm{T}}$ is the control vector of the mobile robot.

The components of the control vector are constrained.

$$u^{-}=-10\le u_i\le 10=u^{+}, i=1,2$$

(35)

A set of initial states is provided.

$${\mathrm{X}}_0=\{x^{0,1},\dots ,x^{0,K}\},$$

(36)

$$x^{0,k}=x^0-{\Delta }_0+{(i)}_3\odot {\Delta }_0, k=1,\dots ,K=24$$

(37)

where $k=j_i-1$, ${(i)}_3$ is a ternary notation of the number $i$ and $\odot$ is the Hadamar product of vectors.

$$j_i=\left\{\begin{array}{l}i, \mathrm{if} i\in \{0,\dots ,11\}\hfill \\ i-3, \mathrm{if} i\in \{15,\dots ,27\}\hfill \end{array}\right., i=0,\dots ,11,15,\dots 26,$$

(38)

$$x^0={[0 0 0]}^{\mathrm{T}},$$

(39)

$${\Delta }_0={\left[5 5 {\displaystyle \frac{\pi }{18}}\right]}^{\mathrm{T}}.$$

(40)

In a sequence of $i$ (38) elements 12, 13 and 14 are deleted because they correspond to the initial conditions $x_1^{0,i}=0$ and $x_2^{0,i}=0$. In total, the initial set includes 24 points of initial states.

The terminal state is specified as

$$x^f={\left[x_1^f x_2^f x_3^f\right]}^{\mathrm{T}}={[0 0 0]}^T.$$

(41)

It is necessary to find a control function (8) that minimizes the following quality criterion:

$$J_3=p_1||x^f-x(t_{f,k},x^{0,k})||+{\displaystyle \sum _{k=1}^{24}{t}_{f,k}}\to \underset{h(x^f-x)}{\min },$$

(42)

where $p_1$ is a penalty coefficient, $p_1=2$,

$$t_{f,k}=\left\{\begin{array}{l}t, \mathrm{if} t<t^{+} \mathrm{and} \left|\right|x^f-x(t,x^{0,k})\left|\right|\le \epsilon \hfill \\ t^{+}. \mathrm{otherwise}\hfill \end{array}\right.,k=1,\dots ,24,$$

(43)

$\epsilon =0.025$, and $t^{+}=1$

$$||x^f-x||=\sqrt{{(x_1^f-x_1)}^2+{(x_2^f-x_2)}^2+{(\sin (x_3^f)-\sin (x_3))}^2}.$$

(44)

According to the above approach, the optimal control problem was solved 24 times, once for each initial state from the set of initial states (36). A direct approach was employed to solve the optimal control problem. The time axis was divided into equal intervals of $\Delta t=0.2$. For each component of control, the parameter values were determined at the boundaries of the intervals, which were then connected by straight lines. As a result, a piecewise linear approximation of the control function was obtained. Considering the constraints on control, the control function had the following form:

$$u_i(t)=v_i(t)=\left\{\begin{array}{l}{u}^{+} \mathrm{if} u^{+}\le {\tilde{u}}_i(t)\hfill \\ u^{-} \mathrm{if} {\tilde{u}}_i(t)\le u^{-}\hfill \\ {\tilde{u}}_i(t), \mathrm{otherwise}\hfill \end{array}\right., i=1,2,$$

(45)

where

$${\tilde{u}}_i(t)=q_{i+(j-1)m}+(q_{i+jm}-q_{i+(j-1)m}){\displaystyle \frac{t-j\Delta t}{\Delta t}},$$

(46)

$$j\Delta t\le t<(j+1)\Delta t, j=1,\dots ,S+1,$$

(47)

$q={[q_1\dots q_{m(S+1)}]}^{\mathrm{T}}$ is a vector of desired parameters, $m$ is a dimension of control, $m=2$, $S$ is a number of time intervals, and

$$S=⌈{\displaystyle \frac{t^{+}}{\Delta t}}⌉=⌈{\displaystyle \frac{1}{0.2}}⌉=5.$$

(48)

The solution to the optimal control problem is the vector of parameters,

$$\tilde{q}(x^{0,j})={\left[{\tilde{q}}_1(x^{0,j})\dots {\tilde{q}}_{12}(x^{0,j})\right]}^T.$$

(49)

To search for the optimal values of the parameter vectors, a hybrid evolutionary algorithm [24,25,26,27] was employed. The obtained parameter vectors are presented in Appendix A.

In the second stage, the control synthesis problem for stabilizing the control object at the terminal point of the state space was addressed using symbolic regression. In solving this problem, the quality criterion specified in criterion (17) was employed to minimize the deviation from the optimal trajectories identified in the previous stage. The control synthesis problem was solved through machine learning utilizing VCBGP. The optimization via the variational genetic algorithm was performed with the following parameter values, where the number of possible solutions in the population was identified as follows: $H=512$, a number of generations $G=128$, a number of possible crossover operation in one generation $R=128$, a depth of variations $d=7$, and a probability of mutation $p_{\mu }=0.75$.

VCBGP obtained the following solution:

$$u_i(t)=h_i(x^f-x)=\left\{\begin{array}{l}{u}^{+} \mathrm{if} u^{+}\le {\tilde{h}}_i(x^f-x)\hfill \\ u^{-} \mathrm{if} {\tilde{h}}_i(x^f-x)\le u^{-}\hfill \\ {\tilde{h}}_i(x^f-x), \mathrm{otherwise}\hfill \end{array}\right., i=1,2,$$

(50)

$${\tilde{h}}_1(x^f-x)=\mathrm{sgn}(A+B)\sqrt{|A+B|}$$

(51)

$${\tilde{h}}_2(x^f-x)=\mathrm{sgn}(\mathrm{sgn}(C)\sqrt{\left|C\right|}+D)\sqrt{|\mathrm{sgn}(C)\sqrt{\left|C\right|}+D|}$$

(52)

$$A=(q_3{({\Delta }_1)}^{-1}{q}_1^2{\Delta }_2\cos ({(\mu ({\Delta }_2)\tanh ({\Delta }_2))}^{3}E),$$

$$B=q_1^3\mathrm{sgn}({\Delta }_3)\sqrt{\left|{\Delta }_3\right|}\exp (\mathrm{sgn}(q_1)\ln (|{\Delta }_3|))+F,$$

$$\begin{gathered} C={(\tanh ({\Delta }_3)({\Delta }_1-{\Delta }_1^3))}^{-1}+{\Delta }_1({\Delta }_2-{\Delta }_2^3)+(-q_2+q_1){({\Delta }_3^3+1)}^3, \\ D=q_1{\Delta }_1+{\rho }_{18}(\sqrt[3]{{\Delta }_1}+\tanh ({\Delta }_3))+G, \\ E={\rho }_{17}({\Delta }_3)\mathrm{sgn}({\Delta }_3)\sqrt{\left|{\Delta }_3\right|}-{\left({\rho }_{17}({\Delta }_3)\mathrm{sgn}({\Delta }_3)\sqrt{\left|{\Delta }_3\right|}\right)}^3, \\ F=q_1{\Delta }_1^2+{\rho }_{17}({\Delta }_3)+{\rho }_{18}({\Delta }_3)-{\left(q_1{\Delta }_1^2+{\rho }_{17}({\Delta }_3)+{\rho }_{18}({\Delta }_3)\right)}^3, \\ G=\exp (-{\rho }_{18}({\Delta }_3){\Delta }_2)\exp (\arctan ({\Delta }_3)+\tanh ({\Delta }_3)), \\ \mu (\alpha )=\left\{\begin{array}{l}\alpha , \mathrm{if} \left|\alpha \right|<1\hfill \\ \mathrm{sgn}(\alpha ), \mathrm{oherwise}\hfill \end{array}\right., {\rho }_{17}(\alpha )=\mathrm{sgn}(\alpha )\ln (|\alpha |+1), \\ {\rho }_{18}(\alpha )=\mathrm{sgn}(\alpha )(\exp (\alpha )-1), {\rho }_{19}(\alpha )=\mathrm{sgn}(\alpha )\exp (-|\alpha |), \end{gathered}$$

${\Delta }_1=x_1^f-x_1$, ${\Delta }_2=x_2^f-x_2$, ${\Delta }_3=x_3^f-x_3$, $q_1=13.65918$, $q_2=10.41113$, $q_3=15.86060$.

Figure 4 and Figure 5 present the codes of the mathematical expressions (51) and (52) in the form of complete binary trees.

In Figure 4 and Figure 5, the numbers in the nodes, except for the leaves of the trees, represent the function numbers with two arguments: 1—$f_{2,1}(z_1,z_2)=z_1+z_2$ and $f_{2,2}(z_1,z_2)=z_1{z}_2$. Numbers near arcs are the function numbers with one argument: 1—$f_{1,1}(z)=z$, 2—$f_{1,2}(z)=z^2$, 3—$f_{1,3}(z)=-z$, 4—$f_{1,4}(z)=\mathrm{sgn}(z)\sqrt{\left|z\right|}$, 5—$f_{1,5}(z)=z^{-1}$, 6—$f_{1,6}(z)=\exp (z)$, 7—$f_{1,7}(z)=\ln (|z|)$, 8—$f_{1,8}(z)=\tanh (z)$, 10—$f_{1,10}(z)=\mathrm{sgn}(z)$, 11—$f_{1,11}(z)=\cos (z)$, 12—$f_{1,12}(z)=\sin (z)$, 13—$f_{1,13}(z)=\arctan (z)$, 14—$f_{1,14}(z)=z^3$, 15—$f_{1,15}(z)=\sqrt[3]{z}$, 16—$f_{1,16}(z)=\mu (z)=\left\{\begin{array}{l}\mathrm{z}, \mathrm{if} \left|z\right|<1\hfill \\ \mathrm{sgn}(z), \mathrm{otherwise}\hfill \end{array}\right.$, 17—$f_{1,17}(z)={\rho }_{17}(z)=\mathrm{sgn}(z)\ln (|z|+1)$ 18—$f_{1,18}(z)={\rho }_{18}(z)=\mathrm{sgn}(z)(\exp (z)-1)$, 19—$f_{1,19}(z)={\rho }_{19}(z)=\mathrm{sgn}(z)\exp (-|z|)$, and 23—$f_{1,23}(z)=z-z^3$. In the leaves, arguments of mathematical equations are presented, ${\Delta }_1=x_1^f-x_1$, ${\Delta }_2=x_2^f-x_2$, ${\Delta }_3=x_3^f-x_3$.

Figure 6 shows robot trajectories from eight initial states (solid black lines), $x^{0,1}$, $x^{0,3}$, $x^{0,7}$, $x^{0,9}$, $x^{0,13}$, $x^{0,15}$ $x^{0,19}$, and $x^{0,21}$, and optimal trajectories (dots) from the same initial states found when solving the optimal control problem.

Figure 7 shows robot trajectories from eight other initial states (solid black lines), $x^{0,4}$, $x^{0,6}$, $x^{0,10}$,$x^{0,12}$, $x^{0,16}$, $x^{0,18}$ $x^{0,22}$, and $x^{0,24}$, and optimal trajectories (dots) from the same initial states. As can be observed, the stabilization system facilitates the motion of the robots in the vicinity of the optimal trajectories.

5. Conclusions

To address the synthesis problem of a stabilization system, machine learning through symbolic regression is employed. To ensure the desired quality properties of the stabilization system, a training sample is utilized during the machine learning process. This training sample is derived from solving the optimal control problem for a set of initial states. The quality criterion in the control synthesis problem is defined as the minimum sum of deviations from all optimal trajectories. This work presents a method for solving the control synthesis problem based on machine learning with supervision. An example of synthesizing a control stabilization system for a mobile wheeled robot with differential drive is demonstrated using symbolic regression techniques. The variational complete binary genetic programming, as a method of symbolic regression, is employed. This symbolic regression method utilizes the principle of small variations in the basic solution and encodes a mathematical expression in the form of complete binary computational trees.

6. Future Work

In future work, it is planned to investigate the influence of stabilization quality on the sensitivity of control systems to external disturbances and inaccuracies in mathematical models. This research aims to enhance the automatic construction of control systems for robots. The investigations will focus on defining the requirements for stabilization systems and determining the values of quality criteria in optimal control problems.