Challenges in control of groups of robots arise from dynamic constraints which assure the absence of collisions between robots. To solve the problem of optimal control for the group of robots we use a two-stage method of synthesis. At the first stage we solve a problem of stabilization of each robot in some point of the state space. At the second stage we search for optimal moving trajectories as sets of points of the state space at which each robot is stabilized. Optimization criteria contain conditions of absence of collisions and other phase constraints. For stabilization we use a symbolic regression method. At the second stage to search for optimal trajectories we use various evolutionary and gradient algorithms of nonlinear programming. In the example, we considered a group of three mobile robots. To search for points of optimal trajectories we used four methods: the genetic algorithm, the particle swarm optimization, the fast gradient descent and the bees algorithm.
The classical approach to solving the applied problem of optimal control for a mobile robot is to solve successively two problems. The first one is a control synthesis problem or ensuring stability to the robot relative to some point of the state space. The second one is the parametric optimal control problem, which consists in finding points of the state space relative to which the control system, synthesized at the first stage, should ensure stability. And sequential switching of the detected points in the state space ensures the movement of the robot from the initial condition to the terminal one, taking into account the phase constraints and the optimal value of the given quality criterion.
The same approach we apply to the control of a group of robots. When solving problems at both stages, we use evolutionary computation methods, for the synthesis problem - the symbolic regression, for the problem of optimal parametric control - evolutionary algorithms.
Problem Statement of Optimal Control of Group of Robots Consider the problem of optimal control for a group of robots.
Given models of the plants
x_ j = f j (xj , uj ), where xj ∈ Rnj , xj = [xj1 . . . xjnj ]T is a state vector of the robot j, uj ∈ Uj ⊆ Rmj , uj = [uj1 . . . ujmj ]T is a control vector of the robot j, Uj is a compact set, j = 1,. . . ,N , N is a number of robots.
Given initial conditions
Given terminal conditions where
xj (0) = x0;j , j = 1, . . . , N.
φjk(xj (tf )) = 0, k = 1, . . . , lj , lj ≤ nj , tf = {t,
if t < t+ and φjk(xj (t)) = 0, k = 1, . . . , lj , j = 1, . . . N ; t+, if t = t+, t+ is a given time limit of control.
Given static phase constraints
αk(xj (t)) ≤ 0, k = 1, . . . , r.
To define dynamic phase constraints, we introduce in the set of robot numbers I = {1, 2, . . . , N } a set of pairs
V = ((i1, j1), (i2, j2), . . . , (iW , jW )), where ik, jk ∈ I, ik ̸= jk.
Then the dynamic phase constraints determining the conditions of closure of pairs of objects are given by βr(xik (t), xjk (t)) ≤ 0, k = 1, . . . , W, r = 1, . . . , R, where W is the number of combinations of N elements taken 2 at a time, W = N !/(2!(N − 1)!) = N (N − 1)/2.
The dynamic phase constraints poses the property of vector commutativity
The points found (10) are stabilization points of robots where the index p increases its value over given intervals of time
Xj = (xj;1, xj;2, . . . , xej;K , εj ), j = 1, . . . , N. e e e
uj = gj (xej;p − xj ), p ← p + (1 − ϑ(tj − t)), tj =← tj + δt(1 − ϑ(tj − t)),
βr(xik (t), xjk (t)) = βr(xjk (t), xik (t)).
J = ∫ tf 0
f0(x1, . . . , xN , u1, . . . , uN ) dt → 0,
It is necessary to find a control function that ensures the movement of robots from the given initial states (2) to the terminal position (3) without violating the phase constraints (5), (7) with the optimal value of the quality criterion (8).
At the first stage we solve the problem of stabilization
uj = gj (xej − xj ), where xj is some point of the state space Rnj .
e To solve the problem, we use one of the methods of symbolic regression [Diveev, 2015a].
At the second stage we find the set of points of the state space and the parameters of switching between them (1) (2) (3) (4) (5) (6) (7) (8) (9) (10) (11) t1 = δt, δt is a given interval of time, ϑ(A) = {1, if A > 0 ;
0, otherwise.
Je = J + ω1h1 + ω2h2,
Coordinates of stabilization points are searched simultaneously for all robots, adding penalties to the target function for violation of phase constraints where ∫ tf ∑r
N ∑ ϑ(αk(xi(t)))dt, h2 = k=1 j=1
0 ∫ tf ∑s
W ∑ ϑ(βk(xj1 (t), xj2 (t)))dt. k=1 w=1 3
Method of Variational Complete Binary Genetic Programming for Synthesis of Stabilization System To solve the problem of control synthesis (9) is finding a multidimensional nonlinear function uj = gj (xej − xj ) that ensures the stability of a system of differential equations
xj = f j (xj , gj (xej − xj )) relative to the point of the state space. We use one of the methods of symbolic regression, the method of variational complete binary genetic programming.
To build a code of binary genetic programming we use the following basic sets: - a set of arguments of mathematical expression - a set of functions with one argument - a set of functions with two arguments
F0 = (q1, . . . , qP , x1, . . . , xN ); F1 = (f1;1(z), f1;2(z), . . . , f1;R(z));
F2 = (f2;1(z1, z2), . . . , f2;S (z1, z2)); - a set of unit elements for functions with two arguments
A set of functions with one argument must include the identity function (12) (13) (14) (15) (16) (17) (18) (19) (21) E2 = (e1, . . . , eM ).
f1;1(z) = z. f2;i(ej , z2) = z2, f2;i(z1, ej ) = z1, Every function with two arguments from (16) has a unit element from the set (17), ∀f2;i(z1, z2) ∈ F2 ∃ej ∈ E2 i ∈ {1, . . . , S}, j ∈ {1, . . . , M }.
To generate a code of binary genetic programming we combine the set of arguments (14) of mathematical expression and the set of unit elements (17) into one ordered set
F = (f1 = q1, . . . , fP = qP , fP +1 = x1, . . . , fP +N = xN , fP +N+1 = e1, . . . , fP +N+M = eM ). (20)
We write down mathematical expression in the form of composition of nested functions and arguments of mathematical expression where f i ∈ F ∪ F1 ∪ F2, i = 1, . . . , K.
Define the inmost nesting depth of an element in the composition. Let the inmost depth is equal L. It means that we have binary computing tree with level L. It has 2L leaves. We write down mathematical expression in the form of a complete binary computing tree (see Fig.1)
In the Fig.1 a is an item number from the set F (20), b is an item number from the set F2 (16), u is an item number from the set F1 (15).
Each level includes numbers of elements from the set F1 and the same numbers of elements from the set F2, and the last level includes numbers from the set F . If the tree has redundant nodes or edges, then we use a number of any function with two arguments for nodes, a number of identity function for edges and a number of unit element for this function with two arguments on the last level. Code of binary variational genetic programming is an ordered set of numbers of elements from the first, second and other levels of the tree C = (u1, b1, u2;1, u2;2, b2;1, b2;2, . . . , uL;1, uL;2, . . . , uL;2L , aL;1, aL;2, . . . , aL;2L ). (22) Consider an example. Let we have the following mathematical expression
y = e qx cos(qx + sin(x)).
For this mathematical expression, we have sets
F1 = (f1;1(z) = z, f1;2(z) = −z, f1;3(z) = ez, f1;4(z) = cos(z), f1;5(z) = sin(z));
F2 = (f2;1(z1, z2) = z1 + z2, f2;2(z1, z2) = z1z2);
F = (f1 = x, f2 = q, f3 = 0, f4 = 1).
Here f3 = 0 is a unit element for addition and f4 = 1 is a unit element for multiplication.
Write down the mathematical expression in the form of composition of functions from these sets y = f1;1(f2;2(f1;3(f1;2(f2;2(f1, f2))), f1;4(f2;1(f2;2(f1, f2), f1;5(f1))))).
The nesting depth in our case is three. The code of the mathematical expression is
C = (1, 2, 3, 4, 2, 1, 1, 2, 1, 5, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 1, 3, 2, 1, 1, 3). Y = (y1, . . . , y6, q, −x, qx, sin(x), q, x, xq, x, q, 0, x, 0, q, x, x, 0, q, 0, x, 0, q, x, x, 0);
Y = (y1, y2, e xq, cos(xq + sin(x)), . . .); Y = (e xq cos(qx + sin(x)), e xq cos(qx + sin(x)), . . .). where w1 is a position in a code, w2 is a new value of the element of the code.
Let is a code of the mathematical expression for level L.
Then and we obtain after the variation (23) a new code
The small variation (23) satisfies the following conditions where i = 1, . . . , L.
A variational genetic algorithm consists of the following stages.
Set a code of some basic solution
Generate a set of ordered sets of small variations where Wi is an ordered set of variations (23) w = [w1 w2]T , C = (c1, . . . , cK ) K = 2L+2 − 2,
w1 w ◦ C = (cz1, . }.|., w{2, . . . , cK ). To search for optimal solution, we use a variational genetic algorithm. It uses principle of small variations of the basic solution [Diveev, 2015b].
Define a small variation of a code of binary analytic programming as an integer vector of two components wi;j = [w1i;j w2i;j ]T , i = 1, . . . , H , j = 1, . . . , l, l is a set of numbers of variations.
We carry out crossover and mutation on the sets of variations. Select two parents Wi = (wi;1, . . . , wi;l) and Wj = (wj;1, . . . , wj;l), and randomly define a point of crossover k ∈ {1, . . . , l} and exchange tails of the parents Wfi = (wi;1, . . . , wi;k 1, wj;k, . . . , wj;l), Wfj = (wj;1, . . . , wj;k 1, wi;k, . . . , wi;l). Define randomly a mutation point µ ∈ {1, . . . , l} and generate new variation in this point µ · wi; w2i; ]T for both new solutions. =
To solve the problem of parametric optimal control, searching for points of the state space Xej = (xj;1, xj;2, . . . , xej;Kj , εj ), we use evolutionary algorithms that work much better than gradient methods for none e unimodal objective function as we show on the experiment. In this paper, we use three popular evolutionary algorithms that are most suitable for solving the parametric optimal control problem: the genetic algorithm [Goldberg, 1989], the particle swarm optimization [Kennedy & Eberhart, 1995]. [Karpenko & Seliverstov, 2010], the bees algorithm [Pham et al., 2006], and the grey wolf algorithm [Mirjalili et al., 2014].
All evolutionary algorithms use a set of possible solutions with a given number H of elements in which the evolution is performed at each iteration a predetermined number W of times to produce new possible solutions. The solution of the problem is considered to be the best possible solution according to the value of the objective function in the resulting set. 5
Consider the task of optimal control of a group of three mobile robots.
Mathematical models of the control objects are the following
x˙ j = uj1 cos θj , y˙j = uj1 sin θj , θ˙j = uj2, where [xj yj θj ]T is a vector of state of the robot j, [uj1 uj2]T is a control vector of the robot j, j = 1, 2, 3.
Given the initial conditions
xj (0) = xj;0, yj (0) = yj;0, θj (0) = θj;0, j = 1, 2, 3.
ui ≤ uij ≤ ui+, i = 1, 2, j = 1, 2, 3. r − √(xj − x )2 + (yj − y )2 ≤ 0, j = 1, 2, 3, where x , y , r are parameters of the phase constraint, r > 0.
Given the terminal conditions
xj (tf ) = xj;f , yj (tf ) = yj;f , θj (tf ) = θj;f , j = 1, 2, 3, tf = t, if t < t+ and 3 ∑√(xj (t) − xj;f )2 + (yj (t) − yj;f )2 + (θj (t) − θj;f )2 ≤ ε; j=1 where
t+, otherwise where t+ and ε are given positive values.
Given the dynamic phase constraints
r˜ − √(x1 − x2)2 + (y1 − y2)2 ≤ 0, r˜ − √(x1 − x3)2 + (y1 − y3)2 ≤ 0, r˜ − √(x2 − x3)2 + (y2 − y3)2 ≤ 0. (35) Taking into account the phase constraints, the quality functional has the following form j=1
j=1 J = tf +
3 ∑ √(xj (t) − xj;f )2 + (yj (t) − yj;f )2 + (θj (t) − θj;f )2 → min .
Je = tf +
3 ∑ √(xj (t) − xj;f )2 + (yj (t) − yj;f )2 + (θj (t) − θj;f )2+ ∫tf 3 0 j=1
∑ ϑ(r − √(xj − x )2 + (yj − y )2) + ϑ(r˜ − √(x1 − x2)2 + (y1 − y2)2)+ +ϑ(r˜ − √(x1 − x3)2 + (y1 − y3)2) + ϑ(r˜ − √(x2 − x3)2 + (y2 − y3)2 dt → min .
In the overall functional (37), the error in fulfilling the terminal conditions and the time period, in which any phase constraints (32), (35) are violated, are added to the time of the control process.
In the computational experiment the parameters of the problem had the following values: x(1)(0) = 10, y(1)(0) = 10, θ(1)(0) = 0, x(2)(0) = 0, y(2)(0) = 10, θ(2)(0) = 0, x(3)(0) = 5, y(3)(0) = 10, θ(3)(0) = 0, u1 = −10, u2 = −10, u1+ = 10, u2+ = 10, x1;f = 0, y1;f = 0, θ1;f = 0, x2;f = 10, y2;f = 0, θ2;f = 0, x3;f = 5, y3;f = 0, θ3;f = 0, x = 5, y = 5, r = 3, r˜ = 2, t+ = 2.8, ε = 0.01, δt = 0.7.
At the first stage, the problem of synthesis of control by the method of variational complete binary genetic programming for one robot was solved. As a result, the following control was obtained ui = ui , if ui < ui ; ui+, if uei > ui+ ; i = 1, 2,
e uei, otherwise, (38) where u1 = sgn(sgn(A)(exp|A| − 1)) ln(| exp |A| − 1| + 1), u2 = B + C, e e A = sgn(ln |q33(q2 − q2 )| + √1 + exp(q3)) × (exp |ln|q33(q2 − q2 )| + √1 + exp(q3))×
3 3 ×sgn(∆x − sgn(∆y)(exp |∆y| − 1)sgn(∆θ)) × ln(∆x − sgn(∆y)(exp |∆y| − 1)sgn(∆θ)| + 1), B = sgn(q3∆θ)(exp |q3∆θ| − 1) + ( µ(∆θ) + 1 − exp(∆y/∆x) ) 1 1 + exp(∆y/∆x) √ C = sgn(q23∆y(∆x + ∆y)) |q23∆y(∆x + ∆y)| × sgn(sgn(∆θ)(exp |∆θ| − 1)(q3 − q33))×
3 ×(ln |(exp |∆θ| − 1)(q3 − q3 )| + 1) + µ(∆x∆y), ∆x = xf − x, ∆y = yf − y, ∆θ = θf − θ, q1 = 0.84180, q2 = 0.65527, q3 = 2.45020.
The obtained control (38) is a nonlinear function uij = gij (xj;1 − xj , yj;1 − yj , θj;1 − θj ), i = 1, 2, j = 1, 2, 3, that depends on the coordinates of the vector of the state of the robot j, j = 1, 2, 3, and ensures the stability of the system of differential equations (29) with respect to the point [xj;1 yj;1 θj;1]T .
At the second stage we solve the problem of parametric optimal control by evolutionary algorithms. It is necessary to find the coordinates of the points [xj;1 yj;1 θj;1]T , [xj;2 yj;2 θj;2]T , j = 1, 2, 3, which are together with the terminal point [xj;f yj;f θj;f ]T the stabilization points of the system (29). Switching from points [xj;1 yj;1 θj;1]T to points [xj;2 yj;2 θj;2]T and from points [xj;2 yj;2 θj;2]T to points [xj;f yj;f θj;f ]T is performed at specified time moments t1 = 1s, t2 = 2s and the points found should provide the minimum value of the functional (37).
When solving the problem of parametric optimal control, the parameters of evolutionary algorithms were chosen so that the total number of calculations of the objective function was approximately the same. The results of the calculations are given in Table 1. For comparison, the problem was also solved by the algorithm of multipoint fast gradient descent, and the best solution of 48 starts was taken. The table shows values of the functional (37) for each of the 10 tests and the number of calculations of the functional in each test. The table contains: GA – genetic algorithm, PSO – particle swarm optimization, BA – bee algorithm, GW – the grey wolf algorithm, MFGD – multipoint fast gradient descent. The last line of the table shows the average values for ten tests.
From the results of the computational experiments it follows that the tested evolutionary algorithms find solutions with approximately the same values of the objective function, and much better than the solutions found by the fast gradient descent. The best solution found is Je = 1.9241: [x1;1 y1;1 θ1;1]T = [0.08539 3.1936 − 0.0469]T , [x2;1 y2;1 θ2;1]T = [5.7851 7.8697 0.5353]T , [x3;1 y3;1 θ3;1]T = [0.0015 1.4369 0.1729]T , [x1;2 y1;2 θ1;2]T = [0.7282 1.5034 − 0.0294]T , [x2;2 y2;2 θ2;2]T = [10.0000 0.5002 0.0927]T , [x3;2 y3;2 θ3;2]T = [9.4238 3.1071 − 0.7933]T . Acknowledgements This work was supported by Russian Foundation for Basic Research, Projects 16-29-04224-ofim, 17-08-01203-a. 2.2117 2.2424 2.4058 2.5942 2.8727 2.5462 2.3855 2.3068 2.1633 2.6650 2.4394