A Genetic Algorithm Incorporating Knowledge and Natural Selection Theory in Adaptive Control of a Wheeled Robot

1. Introduction

Adaptive control methods are still in the mainstream of research in the control of non-linear objects [1,2]. This is due to the fact that control should take into account variable operating conditions. One approach that takes into account variable operating conditions and compensates for the non-linearity of the control object is feedback linearisation [3]. It is based on a model of the dynamics of the controlled object. Therefore, the quality of control obtained in this way strongly depends on the quality of the model, which is affected by the accuracy of identification of its parameters [4,5,6]. It is reduced by the disturbances acting on the object, which can be divided into those that result from modelling inaccuracies and those that are the result of varying operating conditions. In the context of mobile wheeled robots, this can be a change in the mass carried or a change in the resistance to motion of wheel-surface pairs. It therefore becomes crucial to introduce a mechanism into the structure of the control algorithm that allows real-time adjustment of the model parameters. This is done by gradient methods, based on the control error or its function [7,8,9]. This approach integrates adaptive control with optimisation problems [10,11]. Therefore, modern control algorithms make extensive use of artificial intelligence methods, which are based on the search for the minimum function of the model parameters. This opens up the possibility of using models known from machine learning methods, including hybrid algorithms that integrate data-driven approaches with knowledge of the controlled object [12].

In 2023, the IEEE Control Systems Society published a report [13], In 2023, the IEEE Control Systems Society published a report [1] that outlines how control systems will evolve by 2030. Systems will respond to increasing levels of automation across sectors. The report identifies five key methodologies emerging in control systems. One of these is data-driven modelling and control, which can be understood as the integration of machine learning and control systems. A particular expression of this integration is the combination of data-driven models with physics-based models in control algorithms, e.g. in adaptive control. The aforementioned report and the authors of many papers, e.g. [14,15] propose the use of hybrid modelling methods in control which are a combination of data-based and physics-based models. The paper [14] provides an overview of hybrid models used in smart manufacturing. There, a physical model was compared with a data-based model and the benefits of hybrid models were pointed out. An analogous review was included in the paper [16], where 81 articles on the application of models integrating knowledge (physics of the controlled object) and data were selected. As before, reference was made to manufacturing issues. The article [16] provides a comprehensive overview of predictive models used in different industrial sectors. It classifies models based on the level of knowledge used to construct them: full knowledge models, zero knowledge models and partial knowledge models. The paper [15] discusses the development and application of hybrid machine learning methods supported by physics models in the context of cyber-physical systems. The authors emphasise that traditional physics-based and data-only modelling approaches have their limitations. Their combination offers new possibilities for mapping complex physical objects.

It was mentioned earlier that adaptive control algorithms include a real-time parameter estimation mechanism. This process can be understood as an optimisation problem, which consists of selecting the values of the parameters of the mathematical model in such a way as to minimise the control error or its functions. Hence, genetic algorithms, which the literature refers to as nature-inspired approximation optimisation algorithms [17] can be used in the parameter estimation problem. Reviews on genetic algorithms. [18,19,20] and in particular the article [21] state that genetic algorithms belong to the group of metaheuristic algorithms. They are strategies for searching the solution space, typically of an optimisation problem, to obtain an approximate solution.

In this paper, a successful attempt is made to apply a genetic algorithm that is based on the theory of natural selection (GABONST) to the parameter estimation of a mathematical model of a wheeled robot. This algorithm was proposed in the article [22] and is an extension of the classical genetic algorithm by introducing Darwin’s idea of natural selection. Moreover, the proposed algorithm has been extended with knowledge-based mechanisms that address, among other things, the value of the tracking error of the wheeled robot. Another type of knowledge used is that the model parameters are interval constant. In this way, a hybrid algorithm was obtained that combines the features of a genetic algorithm and process knowledge.

In summary, there are three areas of contribution of this paper:

• Synthesis of an adaptive control law that uses a hybrid genetic algorithm to estimate the parameters of a mathematical model of a wheeled robot.
• The use of the GABONST algorithm as a hybrid model based on process knowledge.
• Numerical verification of the proposed adaptive control algorithm for a mobile wheeled robot under varying operating conditions.

The rest of the article, in Section 2, presents a description of the kinematics and dynamics of a mobile wheeled robot. Section 3 presents the synthesis of optimal adaptive control of the robot. The control synthesis was carried out using a mathematical model of the control object. Section 4 presents the genetic algorithm used in the control to optimize the model parameters. Selection, crossover and mutation procedures are presented in detail, as well as how to incorporate physical knowledge of the robot into the genetic algorithm procedures. Section 5 presents and discusses the results of a numerical test, confirming the correctness of the developed theoretical solutions Conclusions on the results of the work are gathered in Section 6.

Translated with DeepL.com (free version).

2. Kinematics and Dynamics of a Wheeled Robot

The mathematical description of robots, in particular wheeled robots, consists of the equations of kinematics and dynamics. The former describe the relationships between the angular velocities of the wheels and the linear velocity of the robot. Their derivation can be based on Denavit-Hartenberg notation, which remains the standard in robotics. The dynamics equations are used to describe and analyse the motion of the robot under the influence of forces and moments acting on it. Their formalism follows from the Lagrange equations and, in the case of mobile robots, the Lagrange equations with multipliers [23,24,25].

2.1. Kinematics

Knowledge of the kinematics equations will allow the angles of rotation and angular velocities of the drive wheels of the mobile robot to be determined. These will be used as the desired trajectory in the tracking task.

Figure 1a shows a scheme of a 2-wheeled mobile robot, consisting of: two driving wheels, denoted by 1 and 2, a self-aligning wheel mounted at point D, and a frame [25]. Its CAD model is shown in Figure 1b. The drive wheels with radius r rotate about their own axis, which does not change position relative to the frame. They are driven by separate drives consisting of electric motors and gears.

The description of the robot’s position and orientation requires three time-varying coordinates, i.e.: $x_H\left(t\right),y_H\left(t\right)$ - the position of the selected characteristic point of the robot (in this case H) in the reference system $x_0{y}_0$, $\beta \left(t\right)$ - the orientation angle of the frame. These coordinates are related by two equations [25]:

$$\left\{\begin{array}{c}{\dot{x}}_H-{\dot{x}}_A-l_3\dot{\beta }sin\beta =0\hfill \\ {\dot{y}}_H-{\dot{y}}_A-l_3\dot{\beta }cos\beta =0\hfill \end{array}\right.$$

(1)

where $v_A$ is the velocity of point A, $l_3$ is the distance between point A and point H and ${\dot{x}}_A\equiv {\dot{x}}_A\left(t\right)\equiv {\displaystyle \frac{dx}{dt}}$, ${\dot{y}}_A\equiv {\dot{y}}_A\left(t\right)\equiv {\displaystyle \frac{dy}{dt}}$, $\dot{\beta }\equiv \dot{\beta }\left(t\right)\equiv {\displaystyle \frac{d\beta }{dt}}$, ${\dot{x}}_H\equiv {\dot{x}}_H\left(t\right)\equiv {\displaystyle \frac{d{x}_H}{dt}}$, ${\dot{y}}_H\equiv {\dot{y}}_H\left(t\right)\equiv {\displaystyle \frac{d{y}_H}{dt}}$.

The projection of the velocity of point A on the axes of the reference system $x_0{y}_0$ is ${\dot{x}}_A=v_{A}cos\beta$ and ${\dot{y}}_A=v_{A}sin\beta$ where $v_A$ is the value of the velocity of point A. Taking this into account, the Equation (1), can be written in the form of

$$\left\{\begin{array}{c}{\dot{x}}_H-v_{A}cos\beta +l_3\dot{\beta }sin\beta =0\hfill \\ {\dot{y}}_H-v_{A}sin\beta -l_3\dot{\beta }cos\beta =0\hfill \end{array}\right.$$

(2)

It is not possible to determine the three coordinates from the two equations, i.e. from (2). Therefore, in order to determine them, a motion path is assumed, which is given by the function

$$f(x_H\left(t\right),y_H\left(t\right))=0$$

(3)

which is differentiable. From the fact that the velocity of a point H is tangent to the path, the equation follows: $grad\overrightarrow{f}\times {\overrightarrow{v}}_H=0$, which can be written in another form, i.e.:

$${\displaystyle \frac{\partial f}{\partial x_H}}{\dot{x}}_H+{\displaystyle \frac{\partial f}{\partial y_H}}{\dot{y}}_H=0$$

(4)

Now the coordinates $x_H$, $y_H$ and $\beta$ can be determined from the following system of equations

$$\left\{\begin{array}{c}{\dot{x}}_H-v_{A}cos\beta +l_3\dot{\beta }sin\beta =0\hfill \\ {\dot{y}}_H-v_{A}sin\beta -l_3\dot{\beta }cos\beta =0\hfill \\ {\displaystyle \frac{\partial f}{\partial x_H}}{\dot{x}}_H+{\displaystyle \frac{\partial f}{\partial y_H}}{\dot{y}}_H=0\hfill \end{array}\right.$$

(5)

For a complete description of the robot’s movement, equations are needed to describe the angular velocities of the drive wheels. These can be written in vector-matrix form

$$\left[\begin{array}{c}{\dot{\alpha }}_1\\ {\dot{\alpha }}_2\end{array}\right]={\displaystyle \frac{1}{r}}\left[\begin{array}{cc}1& l_1\\ 1& -l_1\end{array}\right]\left[\begin{array}{c}{v}_A\\ \dot{\beta }\end{array}\right]$$

(6)

where $l_1=\left|AB\right|=\left|AC\right|$ and ${\dot{\alpha }}_1\equiv {\dot{\alpha }}_1\left(t\right)\equiv {\displaystyle \frac{d{\alpha }_1}{dt}}$, ${\dot{\alpha }}_2\equiv {\dot{\alpha }}_2\left(t\right)\equiv {\displaystyle \frac{d{\alpha }_2}{dt}}$ are the angular velocities of wheels 1 and 2.

By combining the systems of Equations (5) and (6) and solving them, it is possible to determine the five coordinates, i.e. $x_H,y_H,\beta ,{\alpha }_1$ and ${\alpha }_2$ assuming the velocity of point A, i.e. $v_A$ and the path of the robot. Hence, it follows that the mobile wheeled robot is a system with two degrees of freedom.

2.2. Dynamics

An advantageous form of expression for the motion dynamics of a mobile wheeled robot due to control and identification issues is the vector-matrix form

$$\mathbf{M}\ddot{\mathbf{q}}+\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{q}}+\mathbf{F}\left(\dot{\mathbf{q}}\right)=\mathbf{u}$$

(7)

where $\mathbf{q}={[{\alpha }_1,{\alpha }_2]}^T$ is the vector of angles of rotation of the mobile robot’s drive wheels, $\mathbf{u}={[M_1,M_2]}^T$ is the vector of driving moments, $\mathbf{M}$ is the inertia matrix of the robot, $\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{q}}$ is the vector of moments resulting from the centrifugal and Corriolis forces and $\mathbf{F}\left(\dot{\mathbf{q}}\right)$ is the motion resistance vector.

The use of Lagrange equations to describe nonholonomic systems does not directly allow a description of the dynamics in the form of Equation (7) due to the multipliers. Therefore, another mathematical formalism is used in the literature to describe nonholonomic systems, such as the Maggie equation. [25,26].

Taking the parameter vector $\mathbf{a}=[a_1,\dots ,a_6]$, matrices $\mathbf{M},\mathbf{C}$ and the vector $\mathbf{F}$ from the Equation (7) can be written as

$$\mathbf{M}=\left[\begin{array}{cc}{a}_1+a_2+a_3& a_1-a_2\\ a_1-a_2& a_1+a_2+a_3\end{array}\right],\mathbf{C}\left(\dot{\mathbf{q}}\right)=\left[\begin{array}{cc}0& 2a_4({\dot{\alpha }}_2-{\dot{\alpha }}_1)\\ -2a_4({\dot{\alpha }}_2-{\dot{\alpha }}_1)& 0\end{array}\right],\mathbf{F}\left(\dot{\mathbf{q}}\right)=\left[\begin{array}{c}{a}_{5}sign\left({\dot{\alpha }}_1\right)\\ a_{6}sign\left({\dot{\alpha }}_2\right)\end{array}\right]$$

(8)

whereby the $sgn$ function allows for the effect of the direction of movement on the resistance to movement to be taken into account. Due to the discontinuity of the $sgn$ function, its approximation with the $tanh$. function was adopted in the numerical tests. Hence, the vector $\mathbf{F}$ took the form:

$$\mathbf{F}\left(\dot{\mathbf{q}}\right)=\left[\begin{array}{c}{a}_{5}tanh\left(c_F{\dot{\alpha }}_1\right)\\ a_{6}tanh\left(c_F{\dot{\alpha }}_2\right)\end{array}\right]$$

(9)

where $c_F$ is a design constant. The dynamics of the wheeled robot in the form of Equation (7) constitutes its mathematical model used in the numerical tests. Furthermore, it is the basis for the synthesis of the adaptive control law, which is described in the next section.

3. Synthesis of Adaptive Control in a Tracking Task

An adaptive control algorithm is designed in such a way that it can modify its properties under varying operating conditions. In the case described here, the adaptive algorithm is used to control a mobile robot in a tracking task. The task of following a desired trajectory for mobile robots can be understood as determining such a solution for the initial conditions $[\mathbf{q}\left(0\right),\dot{\mathbf{q}}\left(0\right)]\in \Omega =\{(\mathbf{q},\dot{\mathbf{q}}):{\mathbf{q}}_d\left(0\right)=\mathbf{q}\left(0\right),{\dot{\mathbf{q}}}_d\left(0\right)=\dot{\mathbf{q}}\left(0\right)\}$, that converges to a manifold defined as: ${\Omega }_d=\{(\mathbf{q},\dot{\mathbf{q}}):\mathbf{q}={\mathbf{q}}_d,\dot{\mathbf{q}}={\dot{\mathbf{q}}}_d\}$, where ${\mathbf{q}}_d$ is the vector of desired coordinates, ${\dot{\mathbf{q}}}_d$ is the vector of desired angular velocity resulting from the solution of the inverse kinematics problem [25].

The equation describing the dynamics of the mobile robot has the form (7). The tracking error in relation to angle of rotation and angular velocity was defined in the following form

$$\mathbf{e}={\mathbf{q}}_d-\mathbf{q},\dot{\mathbf{e}}={\dot{\mathbf{q}}}_d-\dot{\mathbf{q}}$$

(10)

In addition, a filtered tracking error was defined as

$$\mathbf{s}=\dot{\mathbf{e}}+\Lambda \mathbf{e}$$

(11)

where $\Lambda$ is a diagonal design matrix with positive elements ${\lambda }_i$. By differentiating the generalised tracking error (11), the expression was obtained

$$\dot{\mathbf{s}}=\ddot{\mathbf{e}}+\Lambda \dot{\mathbf{e}}={\ddot{\mathbf{q}}}_d-\ddot{\mathbf{q}}+\Lambda \left({\mathbf{q}}_d-\mathbf{q}\right)$$

(12)

Multiplying both sides of Equation (12) by the inertia matrix $\mathbf{M}$ we obtain

$$\mathbf{M}\dot{\mathbf{s}}=\mathbf{M}{\ddot{\mathbf{q}}}_d-\mathbf{M}\ddot{\mathbf{q}}+\mathbf{M}\Lambda \left({\dot{\mathbf{q}}}_d-\dot{\mathbf{q}}\right)$$

(13)

By substituting Equation (7) into the above equation, the following was obtained

$$\mathbf{M}\dot{\mathbf{s}}=\mathbf{M}{\ddot{\mathbf{q}}}_d+\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{q}}+\mathbf{F}\left(\dot{\mathbf{q}}\right)-\mathbf{u}+\mathbf{M}\Lambda \left({\dot{\mathbf{q}}}_d-\dot{\mathbf{q}}\right)$$

(14)

and after transformation, the result is

$$\mathbf{M}\dot{\mathbf{s}}=\mathbf{M}\left({\ddot{\mathbf{q}}}_d+\Lambda \dot{\mathbf{e}}\right)+\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{q}}+\mathbf{F}\left(\dot{\mathbf{q}}\right)-\mathbf{u}$$

(15)

Multiplying both sides of relation (11) by matrix $-\mathbf{C}\left(\dot{\mathbf{q}}\right)$ yields

$$-\mathbf{C}\mathbf{s}=-\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{e}}-\mathbf{C}\left(\dot{\mathbf{q}}\right)\Lambda \mathbf{e}$$

(16)

Hence it is easy to determine the expression

$$\mathbf{C}\left(\dot{\mathbf{q}}\right)\dot{\mathbf{q}}=-\mathbf{C}\left(\dot{\mathbf{q}}\right)\mathbf{s}+\mathbf{C}\left(\dot{\mathbf{q}}\right)\left({\dot{\mathbf{q}}}_d-\Lambda \mathbf{e}\right)$$

(17)

Taking this into account, relation (15) can be written as a function of the filtered tracking error $\mathbf{s}$

$$\mathbf{M}\dot{\mathbf{s}}=-\mathbf{C}\left(\dot{\mathbf{q}}\right)\mathbf{s}+\mathbf{f}\left(\mathbf{x}\right)-\mathbf{u}$$

(18)

with non-linear function

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{M}\left({\ddot{\mathbf{q}}}_d+\Lambda \dot{\mathbf{e}}\right)+\mathbf{C}\left(\dot{\mathbf{q}}\right)\left({\dot{\mathbf{q}}}_d-\Lambda \mathbf{e}\right)+\mathbf{F}\left(\dot{\mathbf{q}}\right)$$

(19)

whereby the vector $\mathbf{x}=[\mathbf{e},\dot{\mathbf{e}},{\dot{\mathbf{q}}}_d,{\ddot{\mathbf{q}}}_d]$ contains the signals on which $\mathbf{f}\left(\mathbf{x}\right)$ depends. By defining the auxiliary signals

$$\mathbf{\nu }={\dot{\mathbf{q}}}_d-\Lambda \mathbf{e},\dot{\mathbf{\nu }}={\ddot{\mathbf{q}}}_d+\Lambda \dot{\mathbf{e}}$$

(20)

the form of the function $\mathbf{f}\left(\mathbf{x}\right)$ can be simplified to the form

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{M}\dot{\mathbf{\nu }}+\mathbf{C}\left(\dot{\mathbf{q}}\right)\mathbf{\nu }+\mathbf{F}\left(\dot{\mathbf{q}}\right)$$

(21)

Note that the function $\mathbf{f}\left(\mathbf{x}\right)$ can be written in linear form with respect to the model parameters $\mathbf{a}$, i.e.

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{Y}(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\mathbf{a}$$

(22)

The matrix $\mathbf{Y}(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})$ contains the functions resulting from the description of the robot’s dynamics and has the form:

$$\mathbf{Y}(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})=\left[\begin{array}{cccccc}{\dot{\nu }}_1+{\dot{\nu }}_2& {\dot{\nu }}_1-{\dot{\nu }}_2& {\dot{\nu }}_1& 2({\dot{\alpha }}_2-{\dot{\alpha }}_1){\nu }_2& F_1\left({\dot{\alpha }}_1\right)& 0\\ {\dot{\nu }}_1+{\dot{\nu }}_2& -{\dot{\nu }}_1+{\dot{\nu }}_2& {\dot{\nu }}_2& -2({\dot{\alpha }}_2-{\dot{\alpha }}_1){\nu }_1& 0& F_2\left({\dot{\alpha }}_2\right)\end{array}\right]$$

(23)

where ${\nu }_1,{\nu }_2$ are the elements of vector $\mathbf{\nu }$, ${\dot{\nu }}_1,{\dot{\nu }}_2$ are the elements of vector $\dot{\mathbf{\nu }}$, $F_1\left({\dot{\alpha }}_1\right),F_2\left({\dot{\alpha }}_2\right)$ are the elements of vector $\mathbf{F}\left(\dot{\mathbf{q}}\right)$.

Determining the value of the function $\mathbf{f}\left(\mathbf{x}\right)$ requires knowledge of the parameter values. This is difficult to implement, especially in the case of varying operating conditions. Hence, adaptive control uses an estimate of the function $\mathbf{f}\left(\mathbf{x}\right)$ i.e. $\widehat{\mathbf{f}}\left(\mathbf{x}\right)$ which is given as:

$$\widehat{\mathbf{f}}\left(\mathbf{x}\right)=\mathbf{Y}(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\widehat{\mathbf{a}}$$

(24)

where $\widehat{\mathbf{a}}$ is an estimate of the model parameter vector. This approach requires an extension of the control algorithm to include the parameter adaptation law, which is usually implemented as:

$$\dot{\widehat{\mathbf{a}}}=\mathbf{P}{\mathbf{Y}}^T(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\mathbf{s}$$

(25)

It should be noted that this approach requires the assumption that the $\mathbf{a}$ parameters are constant over the time intervals.

Usually, the law of parameter adaptation is applied in continuous time, but in the context of using adaptation techniques operating in discrete time, the law of parameter adaptation given in the form (25) can be written in discrete form [27]:

$${\displaystyle \frac{{\widehat{\mathbf{a}}}_k-{\widehat{\mathbf{a}}}_{k-1}}{h}}=\mathbf{P}{\left[{\mathbf{Y}}^T(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\mathbf{s}\right]}^k$$

(26)

where h is a discrete step. In the context of the application of optimisation methods, the Equation (26) is more favourably written in the form of

$${\widehat{\mathbf{a}}}_k={\widehat{\mathbf{a}}}_{k-1}+\mathbf{P}h{\left[{\mathbf{Y}}^T(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\mathbf{s}\right]}^k$$

(27)

On this basis, the function to be minimised was defined in quadratic form

$$J=argmi{n}_{\widehat{\mathbf{a}}}{\mathbf{E}}^T\mathbf{QE}$$

(28)

where $\mathbf{Q}={\mathbf{Q}}^T>0$ is a design matrix and $\mathbf{E}={\widehat{\mathbf{a}}}_{k-1}+\mathbf{P}h{\left[{\mathbf{Y}}^T(\dot{\mathbf{q}},\mathbf{\nu },\dot{\mathbf{\nu }})\mathbf{s}\right]}^k-{\widehat{\mathbf{a}}}_k$. If optimisation methods, like genetic algorithms, are designed to search for a maximum, the function J can be transformed to the form

$$J^{*}=argma{x}_{\widehat{\mathbf{a}}}{\displaystyle \frac{c}{J+c}}$$

(29)

where $c>0$ is a design constant.

We want to ensure the stability of the closed-loop system without changing the designed compensating control signal $\widehat{\mathbf{f}}\left(\mathbf{x}\right)$ given by Equation (24). This means that the control law $\mathbf{u}$ must be designed in such a way that the filtered tracking error $\mathbf{s}$ is bounded i.e.

$$\left|s_i\right|\le {\Theta }_i,i=1,2$$

(30)

where ${\Theta }_i>0$ is a constant value. For this purpose, the control signal was enriched with a supervising control ${\mathbf{u}}_s$, generated by an additional supervisory term, which is different from zero when the tracking error reaches the limits of the following set:

$$\left\{s_i:\left|s_i\right|\le {\Theta }_i,i=1,2\right\}$$

(31)

In addition, a component responsible for stabilising tracking errors in the form of proportional-differential control expressed as $\mathbf{Ks}$ is included. Finally, to ensure the stability of the closed system, the control law will take the form

$$\mathbf{u}=\widehat{\mathbf{f}}\left(\mathbf{x}\right)+\mathbf{Ks}+{\mathbf{I}}_s{\mathbf{u}}_s$$

(32)

where ${\mathbf{I}}_s=diag[I_{s_1},I_{s_2}]$ is a matrix with elements $I_{s_i}$, which take the value 1 when $s_i⩾{\Theta }_i$ and 0 otherwise. By substituting the Equation (32) into the Equation (18) the result was obtained:

$$\mathbf{M}\dot{\mathbf{s}}=-\mathbf{C}\left(\dot{\mathbf{q}}\right)\mathbf{s}+\mathbf{f}\left(\mathbf{x}\right)-\widehat{\mathbf{f}}\left(\mathbf{x}\right)-\mathbf{Ks}-{\mathbf{I}}_s{\mathbf{u}}_s$$

(33)

By defining a positive definite function $L={\displaystyle \frac{1}{2}}{\mathbf{s}}^T\mathbf{M}\mathbf{s}$, differentiating it and considering the worst-case when tracking errors are large and supervisory control is active, i.e. ${\mathbf{I}}_s=diag[1,1]$, the result is:

$$\dot{L}={\mathbf{s}}^T\mathbf{M}\dot{\mathbf{s}}+{\displaystyle \frac{1}{2}}{\mathbf{s}}^T\dot{\mathbf{M}}\mathbf{s}=-{\mathbf{s}}^T\mathbf{C}\left(\dot{\mathbf{q}}\right)\mathbf{s}+{\mathbf{s}}^T\mathbf{f}\left(\mathbf{x}\right)-{\mathbf{s}}^T\widehat{\mathbf{f}}\left(\mathbf{x}\right)-{\mathbf{s}}^T\mathbf{Ks}-{\mathbf{s}}^T{\mathbf{u}}_s+{\displaystyle \frac{1}{2}}{\mathbf{s}}^T\dot{\mathbf{M}}\mathbf{s}$$

(34)

After transformations, the following form was obtained:

$$\dot{L}=-{\mathbf{s}}^T\mathbf{Ks}+{\displaystyle \frac{1}{2}}{\mathbf{s}}^T\left(\dot{\mathbf{M}}-2\mathbf{C}\right)\mathbf{s}+{\mathbf{s}}^T\left(\mathbf{f}\left(\mathbf{x}\right)-\widehat{\mathbf{f}}\left(\mathbf{x}\right)\right)-{\mathbf{s}}^T{\mathbf{u}}_s$$

(35)

Since $\dot{\mathbf{M}}-2\mathbf{C}$ is a skew-symmetric matrix, the derivative of the scalar function (35) takes the form

$$\dot{L}\le -{\mathbf{s}}^T\mathbf{Ks}+\sum _{i=1}^2\left|s_i\right|\left[|f_i\left(\mathbf{x}\right)|-|{\widehat{f}}_i\left(\mathbf{x}\right)|\right]-\sum _{i=1}^2{s}_i{u}_{si}$$

(36)

If we select the supervisory control signal in the following form

$${\mathbf{u}}_{si}=\left[F_i+\left|{\widehat{f}}_i\left(\mathbf{x}\right)\right|+{\eta }_i\right]sgn\left(s_i\right)$$

(37)

and assuming $|f_i\left(\mathbf{x}\right)|\le F_i,i=1,2,F_i>0,{\eta }_i>0$, we obtain:

$$\dot{L}<-{\mathbf{s}}^T\mathbf{Ks}-\sum _{i=1}^2{\eta }_i\left|s_i\right|$$

(38)

The derivative of the Lyapunov function is negatively definite, so the closed system (18) is stable, for tracking errors $s_i$ such that $|s_i|\ge {\Theta }_i$. When $|s_i|<{\Theta }_i$ then the supervisory control is inactive and the derivative of the Lyapunov function has the form $\dot{L}<-{\mathbf{s}}^T\mathbf{Ks}$. This proves that the filtered tracking error is bounded. The elements of the vector $\widehat{\mathbf{a}}$ are generated by the genetic algorithm, which is formulated in such a way that it searches for solutions in a limited region of the parameter space. Therefore, the estimate of the parameter vector $\widehat{\mathbf{a}}$ is by assumption bounded. Thus, no proof of this fact is required in the stability analysis.

There is a step function $I_{si}$ in the control law (32), which activates the supervisory control signal ${\mathbf{u}}_s$ at the boundary $s_i={\Theta }_i$. To eliminate the high switching frequency of the actuator, an approximation of $sgn\left(s_i\right)$ in the region defined by ${\Theta }_i$ is introduced, according to the relation

$$sgn\left(s_i\right)\approx \left\{\begin{array}{cc}{\displaystyle \frac{s_i}{{\Theta }_i}},\hfill & |s_i|<{\Theta }_i\hfill \\ sgn\left(s_i\right),\hfill & |s_i|\ge {\Theta }_i\hfill \end{array}\right.,i=1,2$$

(39)

A general diagram of the control system is shown graphically in Figure 2.

Finally, the adaptive control of the trailing motion of a wheeled mobile robot is in the form of the Equation (32) and is realised based on the estimation of the non-linear function $\widehat{\mathbf{f}}\left(\mathbf{x}\right)$ and the parameter adaptation law (26), which is implemented as an optimisation problem of a discrete function given in the form (29). In this case, the optimisation algorithm is implemented based on genetic algorithms, which are described in the next section.

4. Classical Genetic Algorithm

A genetic algorithm is a non-linear, discrete and stochastic process, rather than a mathematically controlled algorithm, which was introduced by John Holland of the University of Michigan in 1975 [28]. Since then, this approach has found a permanent place in machine learning methods. Its task is understood as the search for the minimum or maximum of a function in a specified solution space [28,29]. The practice has been introduced that genetic algorithms search for the maximum of the objective function. Furthermore, the paper [28] states that, depending on the particular problem, it may be sufficient to determine a local maximum or a solution that is close to the local or global maximum. Hence the term that the genetic algorithm is an approximation algorithm, as it returns a solution for which the value of the objective function is as high as possible. The search for the maximum of the objective function, which is called the fitness function, is implemented based on the cyclic execution of the selection, crossover and mutation mechanism, which is preceded by the initialisation of the population. From this perspective, the genetic algorithm can be regarded as iterative.

The initialisation of the algorithm, i.e. the generation of the first population in the specified range, appears to have a significant impact on the success of the algorithm. This was confirmed in the paper [19], where the selection of the population was identified as the main challenge determining the quality of the solutions. It is not only its proper selection that is at the heart of initialisation. The choice of the type of chromosome representation is also made at this stage. The papers [20,21] note that the common coding scheme is binary coding, but octal, hexadecimal or decimal coding are also used. The latter was used in this work.

4.1. Selection

Selection reproduces the phenomenon of natural selection, the emergence of a group of the best individuals from a given population, i.e. those for which the values of the fitness function are greatest. The selection operation determines which of the genes in the current generation can constitute the parental group of the next generation. Such an operation is intended to ensure that the optimal solution is approached. Several selection methods are described in the paper [30], including roulette wheel, linear ranked selection or tournament selection.

4.2. Crossover

Individuals that pass the selection process can generate new offspring. Through crossover operators, new offspring are created on the basis of the two parents [31]. Through crossover, using biological comparison, new genetic material is introduced, thereby maintaining genetic diversity [32]. In the context of an optimisation task, the crossover operator prevents the genetic algorithm from converging too quickly and allows a larger region of the solution space to be searched. Many crossover techniques can be found in the literature that differ in the way genes are exchanged between parents. An overview of these methods can be found in the paper [33].

In the work described in this article, the Simulated Binary Crossover (SBX) method is used for crossover. Its operation is based on binary individuals (encoded in decimal), which can be written as vectors $X^{\left(t\right)}=(x_1^{\left(t\right)},x_2^{\left(t\right)},\dots ,x_n^{\left(t\right)})$ and $Y^{\left(t\right)}=(y_1^{\left(t\right)},y_2^{\left(t\right)},\dots ,y_n^{\left(t\right)})$ where t encodes the iteration number of the genetic algorithm. In the SBX method, offspring generation starts with the determination of the coefficient $\beta$, który kontroluje odległość między rodzicami a potomstwem. which controls the distance between parents and offspring. Its value is calculated according to the relation:

$${\beta }_{SBX}=\left\{\begin{array}{cc}{\left(2u\right)}^{{\displaystyle \frac{1}{\eta +1}}},\hfill & u<=0.5,\hfill \\ {\left({\displaystyle \frac{1}{2(1-u)}}\right)}^{{\displaystyle \frac{1}{\eta +1}}},\hfill & \mathrm{if}u>0.5.\hfill \end{array}\right.$$

(40)

where $u\in [0,1]$ is the random crossover coefficient, and $\eta >0$ is the parameter controlling its intensity. Based on the value of the ${\beta }_{SBX}$ coefficient, the values of the offspring genes are determined:

$$\begin{array}{cc}\hfill x_i^{(t+1)}& ={\displaystyle \frac{1}{2}}\left(x_i^{\left(t\right)}+y_i^{\left(t\right)}\right)+{\displaystyle \frac{1}{2}}{\beta }_{SBX}\left(x_i^{\left(t\right)}-y_i^{\left(t\right)}\right),\hfill \\ \hfill y_i^{(t+1)}& ={\displaystyle \frac{1}{2}}\left(x_i^{\left(t\right)}+y_i^{\left(t\right)}\right)+{\displaystyle \frac{1}{2}}{\beta }_{SBX}\left(y_i^{\left(t\right)}-x_i^{\left(t\right)}\right).\hfill \end{array}$$

(41)

Based on the relationship (41) and taking into account the crossover interval (determines which elements of the vector representing the individual will be changed) and the mutation probability, the SBX algorithm can be presented in pseudocode form:

Algorithm 1 Simulated Binary Crossover (SBX)

1: function SBX($X,Y,{\lambda }_{lb},{\lambda }_{ub},p_c\eta$)        ▹ where X,Y is a pair of parents, $\eta$ is the intensity parameter, ${\lambda }_{lb},{\lambda }_{ub}$ is the crossover interval, $p_c$ is the crossover probability 2:     if $rand<p_c$ then 3:         for $i=1$ to n do 4:            if $i\in [{\lambda }_{lb},{\lambda }_{ub}]$ then 5:                choose randomly $u\in [0,1]$ 6:                if $u<=0.5$ then 7:                    ${\beta }_{SBX}={\left(2u\right)}^{{\displaystyle \frac{1}{\eta +1}}}$ 8:                else 9:                    ${\beta }_{SBX}={\left({\displaystyle \frac{1}{2(1-u)}}\right)}^{{\displaystyle \frac{1}{\eta +1}}}$ 10:                end if 11:                $x_i^{(t+1)}={\displaystyle \frac{1}{2}}\left(x_i^{\left(t\right)}+y_i^{\left(t\right)}\right)+{\displaystyle \frac{1}{2}}{\beta }_{SBX}\left(x_i^{\left(t\right)}-y_i^{\left(t\right)}\right)$ 12:                $y_i^{(t+1)}={\displaystyle \frac{1}{2}}\left(x_i^{\left(t\right)}+y_i^{\left(t\right)}\right)+{\displaystyle \frac{1}{2}}{\beta }_{SBX}\left(y_i^{\left(t\right)}-x_i^{\left(t\right)}\right)$ 13:            else 14:                $x_i^{(t+1)}=x_i^{\left(t\right)}$ 15:                $y_i^{(t+1)}=y_i^{\left(t\right)}$ 16:            end if 17:         end for 18:     else 19:         $X^{(t+1)}=X^{\left(t\right)}$ 20:         $Y^{(t+1)}=Y^{\left(t\right)}$ 21:     end if 22: end function

Figure 3 shows the values of the ${\beta }_{SBX}$ coefficient for different values of the $\eta$ parameter.

The ${\beta }_{SBX}$ values for $\eta =100$ are close to 1 which means that the offspring genes will be more similar to those of the parents. In the case of $\eta =1$, the values of ${\beta }_{SBX}$ are more diverse, allowing the generation of more diverse offspring.

4.3. Mutation

The author of dissertation [34] notes, that one of the limitations of genetic algorithms is premature convergence, which results in a narrowing of the evolutionary process to a local optimum. The author of paper [17] sees the cause of this phenomenon in the lack of genetic diversity in the population, which is due to the fact that the crossover operator does not introduce new genes into the population. So this is a problem that needs to be solved. To do this, mutation operators are used, which by changing the values stored in the vectors representing the individuals, new areas of the considered search space can be explored. [34]. A number of mutation techniques can be found in the literature, which differ in the way they work. A list of these can be found in the paper [35]. A large number of citations relating to mutations of individuals encoded by real numbers can be found in the paper [20].

This work uses the Gaussian mutation method, which is one of the most commonly used mutation operators in genetic algorithms with coding by real numbers. It consists in adding a random number generated from a normal distribution to the value of the gene, which allows modifying the values of the elements of the vector $X=(x_1,x_2,\dots ,x_n)$. These are modified according to the equation:

$$x_i^{*}=x_i+\mathcal{N}(0,{\sigma }_t^2)$$

(42)

where $x_i^{*}$ is the value of the gene i after mutation and $\mathcal{N}(0,{\sigma }^2)$ denotes a random number from a normal distribution with mean $\mu =0$ and variance ${\sigma }_t^2$ which can vary depending on the iteration of the genetic algorithm.

The pseudocode of the Gaussian mutation method, which takes into account the mutation probability $p_m$ and the parameter ${\sigma }_t$, can be presented in the form:

Algorithm 2 Gaussian Mutation

1: function GaussianMutation($X,p_m,{\sigma }_t$)        ▹ where X is a chromosome, $p_m$ is the mutation probability, $\sigma$ is the standard deviation 2:     for $i=1$ to n do 3:         choose randomly $u\in [0,1]$ 4:         if $u<p_m$ then 5:            $x_i^{*}=x_i+\mathcal{N}(0,{\sigma }_t^2)$ 6:         end if 7:     end for 8: end function

4.4. Algorithm GABONST

The paper published in 2020 [22] proposes modifications to the classical genetic algorithm by adding mechanisms to incorporate Darwin’s postulates of natural selection. The authors postulate that the genetic algorithm should incorporate a mechanism of natural selection that makes the best individuals likely to survive. These chances for weaker individuals are less, but possible in two cases: 1) crossover with an individual with a high ability to survive in the environment, which can lead to the transmission of favourable genetic features; 2) mutation, which can lead to favourable genetic features. If the individual resulting from these two cases does not meet the requirements of the environment, it may become extinct over time.

This bilogical description of the natural selection method was converted by the authors into the GABONST algorithm, available in the paper [22]. A modified version of it, which was used in this work, is presented below:

• Generate an initial population ${\mathbf{P}}^{\left(0\right)}=\left\{{\mathbf{p}}_1^{\left(0\right)},{\mathbf{p}}_2^{\left(0\right)},\dots ,{\mathbf{p}}_n^{\left(0\right)}\right\}$ consisting of n individuals that are stored as a vector. Set $gen=0$.
• Calculate the value of the fitness function $f_{GA}\left({\mathbf{p}}_i^{gen}\right)$ for each individual in the population ${\mathbf{P}}^{\left(gen\right)}$.
• Return the smallest value of the fitness function for $p_{best}$ percent of the best individuals and store it in the variable $f_{best}$.
• For $i=1,2,\dots ,n$ compare the value of the fitness function $f_{GA}\left({\mathbf{p}}_i^{gen}\right)$ with the value of $f_{best}$:

  (a)

  If $f_{GA}\left({\mathbf{p}}_i^{gen}\right)\ge f_{best}$: mutate an individual using the Gaussian method (with a variance ${\sigma }_1$ and a mutation probability $p_{m_1}$); record in the population ${\mathbf{P}}^{(gen+1)}$;

  (b)

  If $f_{GA}\left({\mathbf{p}}_i^{gen}\right)<f_{best}$:

  • Select a random individual from the set of best individuals $k_{best}$ and write to the variable $\mathbf{best}$.
  • Perform a crossover of an individual ${\mathbf{p}}_i^{gen}$ with individual $\mathbf{best}$ using the SBX method (with parameters ${\eta }_2,p_{c_2}$ for all elements of the individual’s coding vector) and save the single offspring to the variable ${\mathbf{p}}_{cross}$.
  • If $f_{GA}\left({\mathbf{p}}_{cross}\right)\ge f_{best}$ then save the individual ${\mathbf{p}}_{cross}$ in the population ${\mathbf{P}}^{(gen+1)}$
  • If $f_{GA}\left({\mathbf{p}}_{cross}\right)<f_{best}$:
    • Perform a mutation of an individual ${\mathbf{p}}_i^{gen}$ using the Gaussian method (with a variance ${\sigma }_2$ and a mutation probability $p_{m_2}$) and write to the variable ${\mathbf{p}}_{mut}$.
    • If $f_{GA}\left({\mathbf{p}}_{mut}\right)\ge f_{best}$ then save the individual ${\mathbf{p}}_{mut}$ in the population ${\mathbf{P}}^{(gen+1)}$.
    • If $f_{GA}\left({\mathbf{p}}_{mut}\right)<f_{best}$ then generate a new individual and save it in the population ${\mathbf{P}}^{(gen+1)}$.

• Increase $gen$ by 1 and return to step 2.

The above genetic algorithm taking into account the mechanism of natural selection, in the context of this work, is used to adapt the parameters of the estimation of the function $\widehat{\mathbf{f}}\left(\mathbf{x}\right)$ in the control law (32). Therefore, there is no specific stopping criterion. In this application, the GABONST algorithm is executed $N_{GA}$ times within one discrete step. After its completion, the individual whose fitness function is the largest is determined. This is the individual that encodes the value of the parameter estimate $\widehat{\mathbf{a}}$.

4.5. Knowledge in the Genetic Algorithm

Knowledge of the fitness function optimised by the genetic algorithm can be used to control the optimisation process. In the context of the use of the GABONST algorithm in adaptive control, this could be the use of information about the filtered tracking error $\mathbf{s}$, the assumption of equal movement resistance at each driving wheel or the assumption of the variability of the model parameters, which are interval constant. In view of the above statement, several assumptions can be made regarding the assumed parameters of the GABONST algorithm:

• The estimates of the values of the parameters related to resistance to movement in the initial population should be the same and within the assumed range.
• The structure of the GABONST algorithm reflects the assumption of interval constancy of parameters.
• The percentage of individuals, i.e. $p_{best}$ determining how they are processed in the GABONST algorithm, should be dependent on the filtered tracking error $\mathbf{s}$.
• The mutation intensity for the case when $f_{GA}\left({\mathbf{p}}_i^{gen}\right)\ge f_{best}$ should increase as the filtered tracking error $\mathbf{s}$ increases.
• The influence of an individual $best$ on the outcome of a crossover with an individual ${\mathbf{p}}_i^{gen}$ should be reduced as the filtered tracking error $\mathbf{s}$ increases.
• In the case when $f_{GA}\left({\mathbf{p}}_i^{gen}\right)<f_{best}$, the intensity of mutation that can lead to the inclusion of an individual in a new population should increase as the filtered tracking error $\mathbf{s}$ increases.

• The above points, which are written down as suggestions, arise from knowledge of the controlled object. The first point relates to the assumption of an interval of variation in movement resistance, which can be derived from knowledge of the surfaces on which the robot moves. The assumption of equality of movement resistance is related to the assumption that both wheels are on the same surface when the robot starts.

The second and third points refer to the assumption that the parameters of the robot model are constant intervals. This, in turn, means that for a given interval there is no need to adjust the parameters according to the idea of natural selection, but only a small mutation is necessary to fine-tune them. This approach will also reduce the noise associated with the operation of the genetic algorithm. Only if there is a change in the parameter values, which happens at the boundary of the interval, should a larger number of individuals undergo crossover and mutation, which in the GABONST algorithm is done in steps related to natural selection. Whether the robot’s parameters are within the new interval can be determined by the value of the filtered tracking error $\mathbf{s}$, as discussed in point 3.

Points 4-6 deal with mutation and crossover intensity. In the case when $f_{GA}\left({\mathbf{p}}_i^{gen}\right)\ge f_{best}$, an individual ${\mathbf{p}}_i^{gen}$ can be assumed not to survive in a given environment. For the application under consideration, the intensity of mutation and crossover should increase as the filtered tracking error $\mathbf{s}$ increases, on which the parameters of the SBX method (crossover) and the Gaussian method (mutation) can be made dependent. This can be realised based on the error norm of $\mathbf{s}$.

5. Numerical Test

This section describes the numerical test that was performed to verify the correct operation of the proposed adaptive control, in which the knowledge-based GABONST algorithm was used to adapt the parameters of the mathematical model of the wheeled robot. The test assumes that the desired trajectory of the robot point H is a fragment of the Lame curve of the form

$$f(x_H\left(t\right),y_H\left(t\right))={\left({\displaystyle \frac{x_H}{a_L}}\right)}^{2L}+{\left({\displaystyle \frac{y_H}{b_L}}\right)}^{2L}-1$$

(43)

where $a_L,b_L>0$, $L\in {\mathbb{Z}}^{+}$ are the design constants. From the Equations (5) and (6) were determined:

$$\left\{\begin{array}{c}{\dot{x}}_H={\displaystyle \frac{v_A{x}_H{\left(\frac{y_H}{b_L}\right)}^{2L}}{x_{H}cos\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}-y_{H}sin\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}}}\hfill \\ {\dot{y}}_H=-{\displaystyle \frac{v_A{y}_H{\left(\frac{x_H}{a_L}\right)}^{2L}}{x_{H}cos\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}-y_{H}sin\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}}}\hfill \\ \dot{\beta }=-{\displaystyle \frac{v_A{y}_{H}cos\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}+v_A{x}_{H}sin\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}}{l_3\left(x_{H}cos\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}-y_{H}sin\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}\right)}}\hfill \\ {\dot{\alpha }}_1={\displaystyle \frac{v_A}{r}}-{\displaystyle \frac{v_A\left(l_1{y}_{H}cos\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}+l_1{x}_{H}sin\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}\right)}{l_{3}r\left(x_{H}cos\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}-y_{H}sin\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}\right)}}\hfill \\ {\dot{\alpha }}_2={\displaystyle \frac{v_A}{r}}+{\displaystyle \frac{v_A\left(l_1{y}_{H}cos\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}+l_1{x}_{H}sin\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}\right)}{l_{3}r\left(x_{H}cos\left(\beta \right){\left(\frac{y_H}{b_L}\right)}^{2L}-y_{H}sin\left(\beta \right){\left(\frac{x_H}{a_L}\right)}^{2L}\right)}}\hfill \end{array}\right.$$

(44)

Using the Matlab/Simulink package, a coordinate vector $[x_H,y_H,\beta ,{\alpha }_1,{\alpha }_2]$ for a wheeled robot moving along a section of the Lame curve was determined based on the system of Equations (44). The length of the path is determined by the numerical test time ($t_k=25\left[s\right]$) and the profile of velocity $v_A$ which assumed that the acceleration and deceleration phases of the robot occurred in the 3rd and 18th second, respectively, at the maximum steady-state velocity $v_A^{*}=0.4[m/s]$. The other parameters associated with the assumed kinematics equations are: $l_1=0.163\left[m\right]$, $l_3=0.2\left[m\right]$, $r=0.0825\left[m\right]$, $a_L=2\left[m\right]$, $b_L=1\left[m\right]$ and $L=5$.

Figure 4a shows the speed course of the drive wheels, where the acceleration and braking phases are visible. Figure 4b shows the trajectory of the wheeled robot in the reference system $x_0{y}_0$. This is a fragment of the Lame curve, which is described by the equation (43).

The presented trajectory of the robot and the obtained kinematic parameters of the driving wheels will be the basis for determining the filtered tracking error given by Equation (11) in the tracking control task, for a wheeled robot described by the dynamic equation of motion (7). The parameters of the robot were assumed as follows: $a_1=0.0691$, $a_2=0.0511$, $a_3=0.0291$, $a_4=0.0014$, $a_5=a_6=4.3001$. The control law was assumed according to the Equation (32) and for its implementation it was assumed that: $\Lambda =2{\mathbf{I}}_{2\times 2}$, $\mathbf{K}=2{\mathbf{I}}_{2\times 2}$, ${\Theta }_1={\Theta }_2=1.1$, $F_1=F_2=6$, ${\eta }_1={\eta }_2=0.3$.

In the numerical test considered, the model parameters that are related to the resistance to motion, i.e. $a_5$ and $a_6$, were adapted, while the others were assumed constant. For their adaptation, the GABONST algorithm was used, which in this work was enriched with object knowledge. As an initial population, 16 individuals were randomly generated, which are represented by the vectors ${\mathbf{p}}_i=[p_{i1},p_{i2}],i=1,2,\dots ,16$, where $p_{i1},p_{i2}$ is the value of the parameters related to movement resistance. The values of $p_{i1},p_{i2}$ are randomly chosen from the range $[3.6001,5.0001]$. Furthermore, due to the assumption that the wheels of the robot are on the same surface, it was assumed for the initial population that $p_{i1}^{\left(0\right)}=p_{i2}^{\left(0\right)}$.

The objective function that the genetic algorithm maximised is given by the Equation (29) with $c=0.7,\mathbf{Q}=12{\mathbf{I}}_{2\times 2},\mathbf{P}=0.01{\mathbf{I}}_{2\times 2}$. The genetic algorithm was implemented according to the procedure that is presented in Section 4.4. In one discrete step of the control algorithm, which was $h=0.01\left[s\right]$, three loops of the genetic algorithm were implemented. The value of $k_{best}$ was 8. The other parameter values of the GABONST algorithm used in the numerical test were dependent on the knowledge of the object, as indicated in section 4.5. Their detailed values are summarised in Table 1.

In order to force the GABONST algorithm to adapt the model parameters, it was assumed in the numerical test that the parameters related to the movement resistance would have different values in the time intervals. Therefore, it was assumed in the robot model that the parameters $a_5$ and $a_6$ would vary according to the relationship:

$$a_i=\left\{\begin{array}{cc}4.3001\hfill & t\in [0,5)\hfill \\ 5.5901\hfill & t\in [5,15)]\hfill \\ 4.9451\hfill & t\in [15,25]\hfill \end{array}\right.,i=5,6$$

(45)

where t is the simulation time.

The control signals are presented in Figure 5. Figure 5a) shows the waveform of the overall control $\mathbf{u}$ consisting of the stabilising control PD, i.e. ${\mathbf{u}}_{Ks}=\mathbf{Ks}$, the adaptive control $\widehat{f}\left(\mathbf{x}\right)$ and the supervisory controlo ${\mathbf{u}}_s$. For the assumed parameter ${\Theta }_i=1.1$ the supervisory control term ${\mathbf{u}}_s$ did not activate because the condition $|s_i|<{\Theta }_i$ was satisfied. The system remained stable during the movement.

During the initial phase of the movement, the control $\widehat{\mathbf{f}}\left(\mathbf{x}\right)$ (rys. Figure 5b), is immediately activated, as the initial values of the parameter evaluations are consistent with the values assumed in the robot dynamics model and the initial population of the GABONST algorithm was drawn near the nominal values of the parameters $a_5$ and $a_6$. Only when the parameter values in the model are increased by 30 % (which happens at time $t=5\left[s\right]$) is the PD control (Figure 5c) activated. This is related to the higher values of the tracking errors $\mathbf{e}$ and $\dot{\mathbf{e}}$, and therefore also the filtered tracking error $\mathbf{s}$ (rys. 5d). This is due to the maladaptation of the values of the parameter estimates, which are represented by the individual with the maximum value of the fitness function after three loops of the GABONST algorithm in a given iteration step. As the adaptation progresses, the compensating control becomes more and more adequate to the changed dynamics of the object, while the PD control signal decreases. The quality of adaptation is also confirmed by the course of the filtered tracking error $\mathbf{s}$, which remains close to zero after the increase due to the change in the value of the parameters $a_5$ and $a_6$.

The above description is also valid for the situation occurring at time $t=15\left[s\right]$ when the values of the parameters $a_5$ and $a_6$ change and reach 115 % of the nominal values. Then, analogously as before, the role of PD control temporarily increases and then decreases (in amplitude) as the adaptive control $\widehat{f}\left(\mathbf{x}\right)$ adapts to the new values of the changed parameters.

The Figure 6a) shows the course of adaptation of the parameter estimates related to motion resistance. As mentioned earlier, the initial population was chosen so that no significant adaptation is required during the initial phase of the numerical test. This changes when for $t=5\left[s\right]$ there is a change in the value of the parameters $a_5$ and $a_6$. Then, as shown in Figure 6a), the values of the parameter estimates ${\widehat{a}}_5$ and ${\widehat{a}}_6$ adapt to the new values. This is related to the operation of the genetic algorithm, which seeks to maximise the fitness function (Figure 6b). An analogous situation occurs when the parameters are changed for $t=15\left[s\right]$

The above comments on the generalised tracking error $\mathbf{s}$ are confirmed by the waveform of the error $\mathbf{e}$, which values remain small (Figure 7a). Similar conclusions arise from the waveforms of the desired and obtained angular velocities of the robot wheels (Figure 7b). These differ significantly only when there are disturbances, i.e. for $t=5\left[s\right]$ and $t=15\left[s\right]$.

6. Conclusions

In this paper, the synthesis of tracking control of a mobile wheeled robot was carried out. A genetic algorithm that takes into account the idea of natural selection was used to adapt the estimates of parameters required for the implementation of compensating control. The method provides an opportunity to take into account knowledge of the object’s operating conditions, which usually change interval-wise, for example, due to a change in the roadway surface, as studied in the article, or a change in the load being carried. In such a situation, an increase in the intensity of crossover and mutation is required, which is allowed by the algorithm used. This results in a rapid adjustment of the model parameter estimates to the new conditions. On the other hand, when the conditions of the robot do not change significantly, that is, the parameters of the robot do not undergo significant changes, the intensity of crossover and mutation in the algorithm is reduced so that there is no oscillation of parameter estimates.

The results of the simulation indicated the correctness of the assumptions made, in particular that the control system remains stable and the tracking errors are limited. Based on the values of tracking errors, the quality of the tracking task execution should be considered good. In addition, the adapted parameter estimates tend towards the values of the model parameters, and the adaptation algorithm shows sensitivity to parameter changes during the robot’s operation. Further work planned will be the study of the impact of selection, crossover and mutation algorithms on the quality of control with a particular focus on the impact of object knowledge on the parameters of these operations. The next stage of the research will be the verification of solutions using the real robot.