Hierarchical Real-Time Optimal Planning of Collision-Free Trajectories of Collaborative Robots

In collaborative robotics the manipulator trajectory has to be planned to avoid collisions, yet in real-time. In this paper we pose the problem as minimization of a quadratic functional among piecewise linear trajectories in the angular (joint) space. The minimization is subjected to novel nonlinear inequality constraints that simplify the original non-penetration constraints to become cheap to evaluate in real time while still preserving collision-avoidance. The very first and most critical step of the computation is to find an initial trajectory that is free of collisions. To that goal we minimize a weighted sum of the violated constraints until they become feasible or a maximal number of steps is reached. Sometimes an incremental growing of the obstacle helps. By incremental growing we mean that we sequentially solve auxiliary subproblems with obstacles growing from ground or falling from top and use as the initial trajectory the one optimized in the previous step. The initial trajectory is then optimized while preserving feasibility at each step. We solve a sequence of simple-bound constrained quadratic programming problems formulated in the dual space of Lagrange multipliers, which are related to the original linearized inequality constraints that are active or close-to-active. Finally, we refine the trajectory parameterization and repeat the optimization, which we refer to as an hierarchical approach, until an overall prescribed time limit, being well below a second, is reached.


Introduction
Collaborative robotics is an essential topic in the robotic research [1]. A human operator is allowed to support the robotic manipulator with some actions in the working environment [2][3][4]. The sensoric system, typically depth cameras [5], of the robotic manipulator announces a presence of an obstacle. The control system has to be updated with a new collision-free trajectory. This turns out to be a nontrivial task especially as the new trajectory has to be available in much less than a second, which we refer to as real-time.
There is a lot of algorithmic approaches, but basically the problem is split into two. First of all we need to find a feasible, i.e., collision-free trajectory. Secondly, the trajectory is optimized with respect to energy, for instance. For B Dalibor Lukáš dalibor.lukas@vsb.cz 1 VŠB-Technical University of Ostrava, 17. listopadu 15, 708 00 Ostrava-Poruba, Czech Republic the former task, which are referred to as planning problems, variants of graph search algorithms are often employed to arrive at a series of way-points [6,7]. On the other hand optimization problems or nonlinear integrators solve the latter task to follow the way-points in an optimal manner [8][9][10]. When a real-time collision-free planning [11] is required the graph or discrete planning algorithms become tricky to use especially in case of a high-dimensional configuration space. The main novelty of the present paper is to propose for such higher-dimensional problems a new fast planning algorithm implemented in the C programming language. It relies on tailored collision indicators. Additionally our algorithm is adjusted with two novel strategies: an incremental growing of obstacles and a hierarchical expansion of trajectory control points.
Prior to trajectory planning the inverse kinematics [12,13] and dynamics [14,15] of the robotic manipulator is treated. The manipulators may have more degrees of freedom (DOF) than necessarry, which can make the collision-free trajectory planning easier. However, one has to additionally treat singu-lar solutions. In case the DOFs are defined in joints the inverse kinematic problem can be seen as finding real-valued roots of a matrix polynomial. Unfortunately, the numerical stability is not guaranteed. When studying such a system, namely the inverse kinematics of the UR3 manipulator, in this paper we develop a geometric insight into the structure of the roots. To our best knowledge these particular derivations have not been published in literature yet.
The rest of the paper is organized as follows: In Section 2 we put our work in the context of existing approaches. In Section 3 we give a novel and computationaly stable geometric derivation of inverse kinematics of the UR3 manipulator. In Section 4 we set the optimal trajectory problem while presenting a novel collision indicators (constraint functions) that are efficient in the sense that the gradient methods can escape from collisions easily and these collision indicators are fast to evaluate at the same time. In Section 5 we present a collision-free planning algorithm relying on the steepest-descent method, a novel incremental growing of obstacles, and a novel hierarchical optimization strategy. We give numerical validations in Section 6 such that for each given starting and goal positions of the effector an optimized collision-free trajectory is found at 273 milliseconds at latest. The paper is concluded in Section 7.

Related Work
The first part of our paper presents a novel geometric insight into inverse kinematic solutions of a 6-DOF robotic manipulator UR3. We arrive at new analytic formulae that are numerically stable. In literature the systems comprising of joint DOFs are most frequently solved numerically via finding roots of a matrix polynomial. Many methods rely on the Lyapunov stability theory [16]. In [17] the authors employ the Levenberg-Marquart method. Avoiding self-collisions is treated in [18]. Parameterized classes of inverse kinematic solutions for redundant robots are searched for in [19,20]. In [21] the inverse kinematics is tailored to a problem of optimal spraying.
The key part of the collision-avoidance system is the optimal trajectory planning algorithm. First of all an initial feasible trajectory has to be found. Among the most frequently used approaches there are variants of graph algorithms such as a D * -algorithm [22] or an L * -algorithm [23], which is a linear complexity counterpart to the more traditional A * -algorithm. Here edges in the graph represent feasible paths between way-points being the graph vertices. These methods are mostly efficient in the 2-dimensional (2D) motions of robotic vehicles. In higher-dimensional configuration spaces such as 6D in our case, the approach is hardly effective in the real-time regime. Another class of methods that is less effective in higher-dimensions rely on rapidly-exploring random trees [24]. On the other hand among the methods that are efficient in higher dimensions there is a class of potential field methods [20,[25][26][27], in which an artificial potential is introduced such that the trajectory fixed at the start is attracted to the goal and repeled from the obstacles at the same time. In terms of optimization methods [28] this is a penalty approach. A drawback is that it does not always guarantee a collision-free trajectory, in which case the repulsive forces have to be strenghtened. Similar approaches rely on modelling the potential field by the finite element method [29] or by the level-set method [30]. Note that at some specific setups one can develop a geometric approach [31,32] to the trajectory planning. The above mentioned methods can be combined with stochastic approaches [33,34] and neural networks [19,35,36].
In case the optimization function can be well approximated by a quadratic function the Newton methods [20,37] become very efficient as they converge with a quadratic rate and they allow for a bigger trajectory changes towards the optimum. The Newton methods are further used successively in the class of sequential convex or quadratic programming algorithms [19,38,39]. The tricky part of the algorithms is the treatment of collisions. In [40,41] the authors develop a class of interior point methods, where the collision constraints are penalized via so-called barrier functions, which suffer from ill-conditioning and also from the fact that the constraints might be slightly violated similarly to the artificial potential methods. In this paper we opt for an active-set approach, where the feasibility of the trajectory is preserved at each iteration. For a survey of the planning algorithms we refer to [42].
To our best knowledge not many papers are devoted to development of efficient collision indicators meaning that their gradients lead to a fast escape from the collision. This is one of the main topics in the present paper. We rely our collision indicators on a fast approximation formulae for the penetration distance and penetration volumes of mutual intersections of cylinders and intersections of cylinders and voxels representing obstacles. A related approach in literature can be found in [22], where the surfaces are covered with point clouds and interactions of their bounding boxes is treated in a hierarchical approach using binary trees. Similarly in [43] the bodies are approximated by series of spheres.
Important task for the trajectory planning is the trajectory parameterization. Having in mind a subsequent optimal path following, smooth trajectories such as Bézier curves [37,44,45] are preferable. Here the principal question is how many control parameters we shall take. The more we take the closer to the optimum we get, however, the more computational time it spends. Therefore, hierarchical approaches are benefitial. In [46] information about obstacles is updated using a multi-resolution wavelet basis. In [47] the large optimality systems to be solved is replaced by a hierarchy of subsystems, where the obstacles are included successively. Similarly, in [48] the successive optimal subproblems alternate among incorporation of single obstacles. In this paper we rely on a hierarchy of trajectory parameterizations. We search for deformations of sub-optimal trajectories from a previous step. Such deformations were presented in [49] without a hierarchical approach.
Though it is not the subject of this paper for the sake of completeness we shall mention some work devoted to a large and important area of optimal path-following. Here an optimal control problem is solved such that a given feasible trajectory is followed in the most efficient way when considering the robot dynamics. One typically aims at a minimum energy consumption and a minimum torque in joints. We refer to the seminal paper [8], where the unconstrained problem is introduced, and to [9], where collisions are taken into account. We refer to [50] for an actual survey.

Planar Two-Arm Manipulator
The planar two-arm manipulator is the very fundamental setup, which is worth recalling here. We denote the lengths of arms by a 1 and a 2 and the respective rotational degrees of freedom by θ := (θ 1 , θ 2 ) as depicted in Fig. 1 (left). The positions of arm end-points are denoted by x (1) and x (2) := x. They are computed by the simple means of the following forward kinematics: x := a 1 cos θ 1 sin θ 1 For a given position x the inverse kinematics computes the angles θ 1 and θ 2 . It relies on the law of cosines, see Provided that x is reachable, i.e., x ≤ a 1 + a 2 , we get the so-called lefty (L) and righty (R) solutions where γ := arg(x 1 + ı x 2 ). Though the inverse kinematics returns 0, 1 (in case x = a 1 + a 2 ), or 2 solutions, we simplify its notation to IK2 : R 2

Planar Three-Arm Manipulator
Consider a planar three-arm manipulator of the arm lengths a 1 , a 2 , and a 3 with the rotational degrees of freedom θ := (θ 1 , θ 2 , θ 3 ), respectively, as depicted in Fig. 1 (right). The kinematics of the arms is denoted by x (1) , x (2) , and x, respectively, and it is computed as follows: x := a 1 cos θ 1 sin θ 1 where θ := θ 1 + θ 2 + θ 3 . Now the inverse kinematics has θ as a parameter, which is actually prescribed as the direction, where the effector should Hence the position and angle of the third arm is given and we arrive back at the inverse kinematics of a twoarm manipulator. As a mapping the inverse kinematics may again return 0 (in case x > a 1 + a 2 + a 3 ), 1 (in case x = a 1 + a 2 + a 3 ), or 2 (in case x < a 1 + a 2 + a 3 ) solutions, which we simplify to IK3 : The evaluation reads as follows:

Universal Robot UR3
Finally, we consider a 5-arm manipulator with 6 degrees of freedom as depicted in Fig. 2. We denote the arm lengths by a i > 0, i = 1, 2, . . . , 5 and the degrees of freedom by θ 1 , θ 2 , . . . , θ 6 , which are rotations at joints. Additionally, there are three axial displacements d 2 , d 3 , d 4 > 0 of the vertical arms and we denote the total displacement by which is assumed to be positive. At each joint ι ∈ {1, 2, . . . , 5} we introduce a local orthogonal coordinate system described by unit vectors i ι , j ι , k ι ∈ R 3 . The kinematics is computed by the following sequence of formulae for the arm end-points and transformations of the local Obviously, the last degree of freedom θ 6 , the rotation of Joint6, see Fig. 2 (right), does not influence the kinematics unless the robot manipulates an axially nonsymmetric object. Nonetheless, throughout this paper we shall not consider θ 6 any longer. The degrees of freedom of our interest are denoted by θ := (θ 1 , θ 2 , . . . , θ 5 ) ∈ R 5 .
As for the inverse kinematics, we proceed in three steps. First of all, we resolve the fifth arm and θ 1 . Given the effector position x and the unit directional vector v (effector orientation), the fifth arm is determined, We project x 4 onto the plane ρ perpendicular to the direction k 2 := k 3 := k 4 := k 5 The projection reads x ρ 4 := x 4 − d k 2 and x ρ 4 ∈ ρ implies the complex-valued equation (rather than the real system) to be solved for η, which is an auxiliary variable and θ 1 , Provided d ≤ | x 4 |, we get the following two branches of solutions: where Secondly, we resolve the fourth arm by observing that its direction i 4 = i 5 is orthogonal to v as well as to the axial displacement direction k ± 2 . The latter implies that i 4 ∈ ρ, which together with v⊥i 4 gives the two additional branches of solutions, recalling i 1 := (0, 0, 1). Hence, we get up to four projections of x 3 onto ρ, where p, q ∈ {+, −}. This determines the angle θ p,q 5 where j p,q 5 . Finally, θ 2 and θ 3 can be determined from the inverse kinematics of the two-arm planar system in ρ, see Fig. 3 (right). We shift the origin to a 1 i 1 and express x ρ, p,q 3 in the coordinate system i 2 := i 1 and j p 2 , We have where r ∈ {L, R}. The remaining angle can be computed in ρ as follows: To conclude the inverse kinematics of UR3 can be represented as the mapping IKUR3 : where v = 1, p, q ∈ {+, −}, r ∈ {L, R}. In general, there could be either of 0, 1, 2, 4, or 8 solutions.

Optimal Collision-Free Trajectory Problem
Given a starting joint position θ start and ending position and orientation of the effector, x stop and v stop , first of all the related ending degrees of freedom are computed by the inverse kinematics θ stop := IKUR3(x stop , v stop ) so that we do not switch between the branches unless necesarry.
The optimal unconstrained trajectory is the straight line In case of collision we shall search for an as short as possible collisionfree curve.
The parameter domain [0, 1] is decomposed into n + 1 equidistant intervals, over which continuous piecewise linear trajectories are considered so that the straight line is perturbed in a normal direction as depicted in Fig. 4.
The normal directions are spanned by the 4-dimensional space span (n 1 , . . . , n 4 ) where n i = 1. The trajectories pass n break-points and they read where the vector p := ( p 1 , . . . , p n ) ∈ R 4n to be optimized describes perturbations of θ from the unconstrained optimum θ line and where ϕ i (t) are the continuous basis functions that are piecewise linear along the discretization t i := i n+1 , i = 0, 1, . . . , n + 1, with the step h := 1 n+1 , i.e., for i = 1, 2, . . . , n, We denote the 4n-dimensional vector space of the trajectories by T n .
Note that it is rather straightforward to change the nonsmooth piecewise linear basis functions to, e.g., cubic splines or to globally smooth basis functions such as Lagrange polynomials. They can be beneficial in terms of preserving a continuous acceleration in joints. We solve the following nonlinearly constrained quadratic programming problem: 5 6 π ] 4n : g(θ(t)( p))≤0 where the bounds ± 5 6 π are construction limits of joints and where g : T n → R m is the constraint function that avoids collisions.

Constraints Avoiding Collisions
The constraint g(θ(t)( p)) ≤ 0 avoids mutual collisions of arms, collisions of arms with a horizontal workspace construction, to which we further refer to as the ground, collisions of arms and joints with a vertical workspace construction, and collisions of arms and joints with voxels representing the obstacle. Namely, referring to Fig. 2, we consider m := 20 constraints g i as follows: • Constraints i = 1, 2, 3 avoid collisions of Arm1 with Arm4, Arm1 with Arm5, and Arm2 with Arm5. • Constraints i = 4, 5, 6, 7 avoid collisions of Arm2, Arm3, Arm4, and Arm5 with the ground Other collisions do not need to be considered due to the dimensions of the universal robot UR3 and the workspace as described in Section 6.
Notice that arms as well as joints are represented by their cylindrical envelops. The joints the collisions of which do not need to be indicated by a collision of connected arms are Joint2, Joint3, and Joint4. They are exactly those related to axial displacements d 2 , d 3 , d 4 that have not played a role in the inverse kinematics, but they do so now.
To our best knowledge there is no analytic procedure to evaluate collisions of θ (t) for all t ∈ [0, 1]. Therefore we discretize the trajectory parameter interval into N equidistant subintervals and consider the constraints at the worst-case discrete parameter as follows: where g i : R 5 → R.

Mutual Collisions of Arms
We prescribe the constraint as an approximate depth of the penetration of the two arms as follows: where S j , S k are the axes (segments) and r j , r k are the radii of the cylindrical arms. Notice that this approximation of the penetration depth is a strong indicator meaning that it always indicates the true penetration but it can also indicate a false penetration, see Fig. 5 (left). On the other hand it is very cheap to evaluate, which is desired in the real-time optimization. An exact indicator would rely on computing the volume of intersection of the cylinders, which would be computationally demanding.

Collisions of Arms with the Ground
The constraint is prescribed as an approximate depth of the penetration of the arm with the plane x 3 = 0 (the ground) as follows: where a 3 , b 3 are the vertical (third) coordinates of the arm end-points and r is the radius. Again, this constraint also gives positive false indications of collisions, see Fig. 5 (right), but it is cheap to evaluate hence proper for the real-time optimization. An exact constraint would evaluate the volume of intersection of the cylinder with the ground, which would be computationally demanding.

Collisions of Arms and Joints with the Quadrants
The constraint approximates the depth of the penetration of an arm or joint with a quadrant. Consider the horizontal projection of the arm/joint axis, S : } with the axis end-points a, b ∈ R 3 , the radius r , and the first quadrant We define the level-set function The constraint avoiding penetration with the right quadrant Q R := (q 1 .q 2 ) =:q +Q reads while the one avoiding penetration with the left quadrant where S := {(−x 1 , x 2 ) ∈ R 2 : (x 1 , x 2 ) ∈ S}. An exact constraint would evaluate the volume of intersection of the cylinder with the quadrants, which would be computationally demanding.

Collisions of Arms and Joints with the Obstacle
This constraint approximates the depth-times-volume of penetration of an arm or joint into the voxels representing the obstacle. The product of depth and volume turned out to resolve more collisions than the single depth constraint, like 23 and 24, or the single volume constraint. The voxels are received from cameras and they typically envelop hands of a human operator. Consider a cylindrical arm/joint S determined by radius r and axial end-points a, b. Further, consider the union of all the voxels. We compute the constraint in the following four steps: 1. Detect voxels boundary := ∂ . It collects those square voxel faces that are not shared with another voxel. This step has a quadratic complexity in terms of numbers of voxels, but it is performed only once per evaluation of g.

Compute intervals of penetrations of the four edges x ±i (t)
into , This is done by computing all intersection points of x ±i (t) with and sorting them along the parameter t. The consecutive pairs of the sorted parameters determine intervals inside, which are of our interest, and outside . 4. Compute the constraint by summing up penetrations of the four edges where |I k j | is the length of I k j and t k j is the mid-point of I k j .

Finding an Initial Collision-Free Trajectory
The first and most crucial step is to find a feasible trajectory. A natural idea would be to search for a feasible polygonal trajectory in the 5-dimensional space of joint degrees of freedom. However, the global search suffers from an exponential computational complexity. Hence we rely on the following local search algorithm: Given a (zero) non-feasible initial trajectory p 0 , we proceed the steepest-descent iterations until p k+1 is feasible or k, and hence computational time, reaches a given maximum. The descent direction is the negative weighted gradient of non-feasible constraints, and α k is computed by the line-search (bisections) method The bisection starts at the maximal feasible If we do not succeed we adopt an incremental growing of the obstacle, represented by the voxels, from the ground as sketched in Fig. 6. First of all we find the maximal vertical coordinate x 3 := arg max x∈ x 3 . Starting at the unconstrained optimum p := 0, we solve a sequence of M problems to find a feasible trajectory by algorithm 31, 32, 33, while incrementaly growing the obstacle as follows: The feasible path of the m-th problem is a starting point to the (m + 1)-st problem. If we still do not succeed, we let the obstacle fall from the top as follows: Fig. 6 Incremental growing of the obstacle from the ground

Optimization
After we find a collision-free trajectory, denoted again by p 0 , g( p 0 ) ≤ 0, we finally continue to optimize it. We solve the nonlinearly constrained quadratic programming problem 21 approximately (numerically) by means of the presented optimization algorithms. We employ the sequential quadratic programming and the active-set approach, the iterations of which read where the following quadratic programming subproblem is solved: The constraint matrix G k A comprises of columns ∇g i ( p k ), i ∈ A, at active or close-to-active indices. These are chosen at the 10% level of the maximal constraint, i.e., Moreover, linearly dependent columns of G k A are removed. Hence, 38 is a well-posed linear saddle-point problem that reads as follows: The latter admits the solution

Hierarchical Optimization
For higher numbers of trajectory break points n, i.e., higher dimensions 4n, finding a feasible design and its optimization become difficult. Therefore we employ the idea of multigrid methods [51] to solve the problem through a nested hierarchy. At the initial level l := 1 we start with n (1) := 1 and the unconstrained optimum p (1) 0 := 0. The resulting, so-called triangular, trajectory p (1) * , no matter whether it is feasible or not, is taken as the initial guess p (2) 0 at the next level l := 2. We introduce two additional break points as depicted in Fig. 7, hence n (2) := 3. The process continues at this level and we arrive at the, so-called pentagonal, trajectory p (2) * . The latter is interpolated to the next level l := 3, n (3) := 7, where we stop. The resulting trajectory is called nonagonal.
For our choice of basis functions 20 we arrive at the following interpolation matrices P (1,2) and P (2,3) : where I is the 4-by-4 identity matrix and ⊗ is the Kronecker product meaning that each entry of the matrix on the left is multiplied with I. The columns (. . . , 1/2, 1, 1/2, . . . ) T are coordinates (shapes) of the basis functions ϕ (l)

Numerical Experiments
First of all we shall determine the sizes of the UR3 manipulator and parameters of the optimization problem 21. In the following all the sizes are in meters and the degrees of freedom are in radians. The parameters of UR3 are as follows: The steepest-descent optimization algorithm relies on forward numerical derivatives with the step 10 −6 (radians). The number of steps of the incremental growing is M := 10.
We illustrate performance and behaviour of the algorithm on a model situation depicted in Fig. 8.
We consider the following joint coordinates (in radians) of the manipulator:  where the latter corresponds to the ending position of the effector The obstacle G is a cube, discretized into 4 3 = 64 voxels, which is placed at the center at c := (−0.075, −0.6, 0.325). The size of the cube edges is 0.2. The cube is situated parallel to the coordinate system.
In Tab. 1 we give results of numerical simulations for the following 27 shifts of the center of G, where the factor 0.1 [m] is a relative shift of the obstacle that was chosen to simulate as many collisions as possible: In all these simulation the algorithm found an optimal trajectory unless there was a collision with the obstacle at the ending time instance t stop := 1, which was indicated in less than a millisecond. In some cases the straight line was free of collisions, which was indicated at 1 millisecond. In the other cases typical total computational time was in tens of milliseconds if there was no change of the inverse kinematics branch and below 300 milliseconds if there was a change of the branch.
The first column of Table 1 indicates the shift of the obstacle center. The second column displays whether the optimization succeeded or not and for which inverse kinematic configuration of θ stop . The configuration 0 denotes the unchanged branch. At some cases the branch had to be changed to the second (nearest to unchanged) one, whichis  Examples of a real collision-avoidance in our in-house digital-twin system denoted by 1. A change to the third branch has never happened. In the third column we depict the number m of growing steps from the ground and from the top. For instance the entry '3,3,3/19,19,19' means that the algorithm failed at step m := 3 of the growing from the ground 35 at all the levels l = 1, 2, 3 as well as it failed at step m := 19 of the growing from the top 36 at all the levels l = 1, 2, 3. Recall that M := 10. If the algorithm found a feasible path without the growing it is indicated by 'no need'. The entry '7,10/-' indicates that the algorithm failed with the growing from the ground at step m := 7 at level l = 1, but then it succeeded at level l = 2 so that the growing from the top was not needed, which is indicated by that '-'. In the fourth column of Table 1 the level at which we found optimum is displayed. The fifth column displays the number of steps needed to find a feasible trajectory. The sixth column displays the number of the forthcoming steps to find an optimal trajectory. The last column displays the total computational time of the simulation. Finally, we shall note that the algorithm presented in this paper is a part of a real system including the UR3 manipulator and cameras. In Fig. 9 we sketch a typical human-avoidance situation performed in this real system. The system is implemented in C++. It receives data from a sensory subsystem over Ethernet via UDP (User Datagram Protocol) in the form of a set of voxel centers representing dynamic obstacles. Then the trajectory is updated and the physical UR3 robot is controlled using the RTDE (Real-Time Data Exchange) protocol (Ethernet-based) provided by Universal Robots. Collaborative robots by other manufacturers could be also used, which would require a different communication protocol. 3D visualization is done using our in-house graphical engine based on DirectX 11.1. The camera system relies on an in-house dedicated control unit. It is described in [52].

Conclusion
In this paper we have presented a novel approach to optimal planning of collision-free trajectories of robotic manipulators. We arrive at a code that is capable for real-time usage. Namely, the trajectory is found in a quarter of second at latest. Our method relies on • novel numerically stable formulae derived by geometrical arguments, • novel fast approximations of collision indicators relying on penetration depth and penetration volume of mutual intersections of cylindrical arms and intersections of cylindrical arms with voxels representing obstacles, • a novel hierarchical approach within the steepest-descent optimization algorithm.
The ongoing research is now devoted to development of a novel hierarchical optimal path-following method.
Author Contributions Dalibor Lukáš has derived the inverse kinematics, proposed the novel collision-free indicators, carried out the mathematical research, and implemented code of the novel optimization methods. Tomáš Kot has provided robot's specification, software framework for the numerical simulations, and the physical environment in a lab.
Funding Open access publishing supported by the National Technical Library in Prague. This work was funded by the Research Platform focused on Industry 4.0 and Robotics in Ostrava Agglomeration project, project number CZ.02.1.01/0.0/0.0/17_049/0008425 within the Operational Programme Research, Development and Education.
Code Availability Not applicable.

Competing Interests
The authors have no relevant financial or nonfinancial interests to disclose.
Ethics Approval Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. In 2020, he finished his habilitation at the same university, where he currently works as a Senior Researcher. His research activities focus on complex simulations and control of mechatronic systems, visualisation, application of virtual and augmented reality in robotics, optimisation of layouts of robotised workplaces, algorithms for automatic design of an optimal kinematic structure of a robotic manipulator suitable for a given task, and lately also collision avoidance for collaborative robots sharing workspace with human workers.