Antun Skuric

Vincent Padois

David Daney

Approximating robot reachable space using convex polytopes

Keywords: Reachability analysis, Convex polytopes, Collaborative robotics

This paper presents an approach for approximating the reachable space of robotic manipulators based on convex polytopes. The proposed approach predicts the reachable space over a given time horizon based on the robot's actuation limits and kinematic constraints. The approach is furthermore extended to integrate the robot's environment, assuming it can be expressed in a form of linear constraints, and to account for the robot's link geometry. The accuracy of the proposed method is evaluated using simulations of robot's nonlinear dynamics and it is compared against the cartesian space limits, usually provided by manufacturers in standard datasheets. The accuracy analysis results show that the proposed method has good performance for the time horizons up to t h ≤250ms, encapsulating most of the simulated robot's reachable space while maintaining comparable volume. For a 7 dof robot, the method has an average execution time of 50ms, independent of the horizon time, potentially enabling real-time applications.

Introduction

Collaborative robots, designed for safe physical interaction with humans, have a great potential to find their way in many spheres of industry, research and potentially become a part of our everyday lives. Apart from replacing rigid and dangerous manipulators in the industrial applications, their main potential is creating new robotic assistance scenarios leveraging high degree of physical human-robot interaction. Such scenarios demand not only many safety guarantees, but a high degree of operator's situation awareness both in terms of robot's behaviour and its physical abilities in real-time.

A metric unifying the robot's physical abilities and its movement capacity, important for safety and performance analysis, is called the robot's reachable space [START_REF] Althoff | Reachability analysis of nonlinear differential-algebraic systems[END_REF]. It represents the space of reachable robot's positions (ex. end-effector cartesian positions) for a certain time horizon, given different assumptions on its physical abilities. As robotic systems are highly non-linear and featuring complex dynamics, the true reachable space of robots is very complex to characterise and calculate. Therefore, approximation techniques are necessary in order to yield practical solutions. Pereira and Athhoff [START_REF] Pereira | Calculating human reachable occupancy for guaranteed collision-free planning[END_REF] have developed an interval analysis based approach, later improved by Shrepp et al. [START_REF] Schepp | Sara: A tool for safe human-robot coexistence and collaboration through reachability analysis[END_REF], to approximate the human arm reachable space, modeled as a serial robotic manipulator, using a set of spheres and cylinders. However, in many cases the approximation of the reachable space using these shapes is impractical due to their non-linear nature.

In this paper, an approximation approach of the robot's reachable space using convex polytopes is proposed, leveraging several of their key characteristics. Convex polytopes can be represented as a set of linear inequalities Ax ≤ b, which may be directly used with different optimisation techniques. Additionally, this set of linear inequalities can be intuitively extended with different environmental and user defined constraints. Since the reachable space is often low dimensional (≤3D) these polytopes can be easily visualised, having the structure of a triangulated mesh, and potentially provide a human operator with valuable information about the robot's capacity. Furthermore, operations over polytopes such as Minkowski sum and intersection are well defined and efficient to calculate, enabling for the intuitive extension of the proposed method to account for robot's link geometry. Finally, the efficiency of polytope enumeration techniques for low dimensional polytopes has a potential to make this metric real-time capable.

Even though convex polytopes are widely used to characterise different metrics of robot's physical abilities [START_REF] Skuric | On-line force capability evaluation based on efficient polytope vertex search[END_REF]: force, acceleration and velocity, to the best of our knowledge they have not yet been used to characterise the reachable space of a robotic manipulator. The proposed approach is partially inspired by the works of Long et al. [START_REF] Long | Evaluating robot manipulability in constrained environments by velocity polytope reduction[END_REF] on the constrained manipulability (velocity) polytopes.

Depending on whether the application requires an under or over approximation, or perhaps time-efficiency, the evaluation metrics of the accuracy of the reachable space approximation will differ significantly. In this paper, in addition to the execution time analysis, three accuracy measures are chosen to enable the evaluation of the proposed approach's performance and study its limitations.

The paper has the following structure. In the chapter 2 the formal definition of the approach is introduced, followed by the chapter 3, describing the procedure of benchmarking and the accuracy analysis. In the chapter 4, the results of the accuracy and execution time analysis of the approach are given and the chapter 5 discusses the main limitation as well as potential applications.

Convex polytope of reachable space

Robot dynamics in joint space can be expressed as M (q) q + C(q, q) q + τ g (q)

τ d = τ (1)
where M is the mass matrix, C is coriolis matrix and τ g is the gravity vector. τ is the applied joint torque vector and τ d is the equivalent joint torque vector due to the coriolis and gravity effects.

For a robot with n actuated degrees of freedom (DOF), the robot's joint torque τ , velocity q and joint angles q are n-dimensional vectors, limited by the robot's hardware

1 τ ∈ [τ min , τ max ] , q ∈ [ qmin , qmax ] , q ∈ [q min , q max ] (2)
For a given moment in time k and for given robot configuration q k , the affine relationship between joint torques τ and accelerations qk can be expressed as

qk = M -1 (q k )(τ -τ d ) = M -1 k (τ -τ d ) (3)
Considering fixed joint acceleration qk during given horizon length t h , an approximation of the robot's joint velocity qk+1 and position q k+1 , at the end of horizon, can be calculated using numerical integration (forward Euler method)

qk+1 = M -1 k t h (τ -τ d ) + qk , q k+1 = M -1 k t 2 h 2 (τ -τ d ) + qk t h + q k (4)
This linear numerical integration (4) considers the robot's dynamics [START_REF] Althoff | Reachability analysis of nonlinear differential-algebraic systems[END_REF], and the applied joint torque τ , fixed during the horizon time t h . Such assumption is reasonable only for short horizon times t h which in term represents a limitation of the proposed method.

The relationship between the m-dimensional task space velocity and acceleration of certain frame on the robot (for example end-effector frame) and the n-dimensional joint space equivalents is defined through the corresponding jacobian matrix J(q) and its time derivative J(q). ẋk

= J(q k ) qk = J k qk , ẍk = J k qk + Jk qk (5)
Finally, given the horizon of interest t h , the predicted cartesian position x k+1 can then be expressed as

x k+1 = ẍk t 2 h 2 + ẋk t h +x k = J k M -1 k t 2 h 2 τ -J k M -1 k t 2 h 2 τ d ∆x k,dyn + Jk qk t 2 h 2 + ẋk t h ∆x k,vel +x k x * k+1 (6)
Where the x * k+1 = ∆x k,dyn +∆x k,vel + x k is a predicted position vector based on the joint configuration q k , joint velocity qk (or cartesian velocity ẋk = J qk ) and the joint torque τ d in the current time-step k, while the first term describes the influence of the applied joint torques τ .

The convex polytope of reachable cartesian space P x can then be defined as a set of all the possible cartesian positions x k+1 at the end of the horizon t h , achieved by any combination of joint torques τ within robot's physical limitations [START_REF] Bowling | The dynamic capability equations: a new tool for analyzing robotic manipulator performance[END_REF], given the robot's current state qk , q k and τ d . ). Third image from the right shows the robot in its initial position q = [0, 0, 0, -π/2, 0, 3π/5, 0], close to the center of the joint ranges. In the second image form the right, the robot's joints q1 and q3 at their limits, q = [0, -1.59, 0, -2.9, 0, 3π/5, 0, 0]. On the last image, robot's joints q0 and q2 are at their limits q = [-2.72, 0, -2.72, -π/2, 0, 13π/5, 0, 0], preventing the robot to rotate around the z axis in one direction. Fig. 2. Three images on the left show the shifting effect on the reachable space polytope Px produced by certain cartesian velocity at the beginning of the horizon, from left to right: ẋk = [0, 0, 0], [0, 1, 0] and [0, 0, -1]m/s. Three images on the right show the reducing effect on the polytope Px by three different carried object masses (left to right) mo = 0, 2 and 5kg. Horizon used is t h =0.15s and the robot is in initial configuration.

P x = {x k+1 ∈ R m | x k+1 = J k M -1 k t 2 h 2 τ + x * k+1 , τ ∈ [τ min , τ max ] , M -1 k t h (τ -τ d ) + qk ∈ [ qmin , qmax ] , M -1 k t 2 h 2 (τ -τ d ) + qk t h + q k ∈ [q min , q max ]} (7)
Figures 1 and2 showcase the influences of different horizon lengths t h , robot constraints and robot movement ( qk , τ d ) on the shape of P x . Chapter 2.4 describes an approach to enumerate the reachable space polytope.

Influence of the carried object

If the robot is carrying a payload, an object with a mass m o and inertia i o attached to its end-effector, this object will have an influence on the robot's dynamics, modifying the mass matrix M , coriolis matrix C and the gravity vector τ g [START_REF] Hamad | The role of robot payload in the safety map framework[END_REF]. The augmented dynamical model of the robot will have reduced acceleration capabilities due to the added effort necessary for the object's movements. The figure 2 shows the influence of different object masses m o , on the resulting reachable space polytope P x for the Franka Emika Panda robot.

Integration of the environment

If the robot's environment can be defined as a set of convex inequalities

A e x ≤ b e (8)
these constraints can be directly transformed to the constraints of the joint torque τ using the equation ( 6)

A e x k+1 = A e J k M -1 k t 2 h 2 (τ -τ d ) + A e x * k+1 ≤ b e (9)
and then included in the P x calculation ( 7)

P x = {x k+1 ∈ R m | x k+1 = J k M -1 k t 2 h 2 τ + x * k+1 , τ ∈ [τ min , τ max ] , M -1 k t h (τ -τ d ) + qk ∈ [ qmin , qmax ] , M -1 k t 2 h 2 (τ -τ d ) + qk t h + q k ∈ [q min , q max ] , A e J k M -1 k t 2 h 2 (τ -τ d ) + A e x * k+1 ≤ b e } (10)
Figure 3 showcases the influence of the environmental constraints on the shape of the P x .

Integration of robot's link geometry

If the robot link l i is considered to be a straight line, to calculate the predicted reachable space of a robot's link l i the polytope P x can be calculated at the start P xs and the end P xe of the line l i . Finally the reachable space of this idealised robot link l i can then be calculated as the convex-hull (CH) of the polytopes P xs and If, instead of a straight line, a space captured by a robot's link l i can be expressed as a convex set of constraints or a polytope

P xe P xli = CH (P xs , P xe ) (11)
L i = x l | A l x l ≤ b l (12)
the reachable space polytope P xli can be extended to account for this link geometry by first finding the vertices of the polytope L, followed by computing the reachable space polytope P x for each one of the k vertices, and finally computing their convex hull.

P xl,i = CH (P x1 , . . . , P xk ) (13)
By evaluating the polytope P xli for each robot's link one can obtain an efficient approximation of the envelope of the reachable space of the entire robot, as shown on figure 4.

Enumerating reachable space polytope

The reachable space polytope P x belongs to the family of over-determined problems

P x = {x ∈ R m | x = P τ , Aτ ≤ b }, τ ∈ R n , n ≥ m (14)
where the inequality constraint Aτ ≤ b encapsulates all the robot constraints (2) as well as the environmental constraints (9), and the matrix

J k M -1 k t 2 h
2 becomes the projection matrix P . Vector τ corresponds to the joint torque vector and the vector x corresponds to the difference x k+1 -x * k+1 . In order to fully define and simplify this formulation, the set of all the vertices (V-representation) or all the half-spaces (H-representation) of this polytope has to be found. Vertex representation of a polytope is a list of all of its vertices v i and is commonly used for visualisation purposes V = {v 0 , v 1 , v 2 , . . . }, while half-space representation is defined as a list of inequalities Hx ≤ d, each corresponding to one of the polytope facets. This matrix inequality equation is more suitable for applications involving different optimisation strategies.

The reachable space polytope ( 14) is a particular case of the projection of a high-dimensional polytope onto a lower dimensional space. Therefore, the set of constraints Aτ ≤ b form a H-representation of a polytope P τ in the ndimensional joint torque space.

P τ = {τ ∈ R n | Aτ ≤ b } ( 15
)
This polytope is then projected using the projection matrix P into the mdimensional cartesian space, usually much lower dimensional, forming the reachable space polytope P x .

P x = {x ∈ R m | x = P τ , τ ∈ P τ } (16)
The most straight forward way of approaching the enumeration of this polytope is to first enumerate the vertices of the polytope P τ , project them to the cartesian space using the projection matrix P , followed by a Convex-Hull algorithm to extract the real vertices v i of the polytope P x from the projected points and its half-space representation.

In order to avoid enumerating all the vertices of the high-dimensional joint torque polytope P τ , as the most complex operation of this polytope enumeration, Iterative convex hull method [START_REF] Skuric | On-line feasible wrench polytope evaluation based on human musculoskeletal models: An iterative convex hull method[END_REF] (ICHM) is used, allowing to directly enumerate low-dimensional cartesian polytope P x . This iterative method is based on successively applying Linear programming (LP) and Convex-Hull (CH) algorithms and is able to find both vertex and half-space representation of this polytope at the same time. ICHM algorithm is defined for a family of sets:

{ x ∈ R m | Ax = By, Cy ≤ d }, y ∈ R n , m ≤ n (17)
Therefore the usage of the ICHM algorithm for enumerating the reachable space polytope is rather straight forward, by setting its matrix A to identity A = I m×m , matrix B becomes the projection matrix P , the inequality constraint Cy ≤ d then becomes the set of the constraints Aτ ≤ b.

Analysing the approximation accuracy

Defining metrics of interest for accuracy analysis is highly dependent of the applications these methods will be used on. If the application requires an underapproximation of the reachable set or an over-approximation, the metrics of interest will not be the same. In this paper, the algorithms accuracy is analysed using three different quantitative metrics, each one providing an insight in different characteristics and limitations of the method.

Additionally as this approximation method has a potential to be used for real-time applications, execution time analysis is performed as well.

In the extent of the proposed experiments, the considered reachable space is the reachable space of the robot's end effector. For each vertex of the Px, joint torque τ generating the vertex is applied to the robot simulation (figures 2 and 3). All the robot's end-effector positions in all the simulation time steps ∆t are saved for further analysis.

Metrics definition

To analyse the approximation accuracy of the reachable space polytope P x , a simple discrete non-linear robot dynamics simulation, subject to the constraints (2) is employed. Simulation is carried out by discretizing the equation ( 1) to determine the instantaneous joint acceleration qk and followed by its numerical integration to determine joint velocity qk+1 and position q k+1 in the step k+1.

qk = M -1 (q k ) (τ -C(q k , qk ) qk + τ g (q k )) qk+1 = qk ∆t + qk q k+1 = qk ∆t 2 2 + qk ∆t + q k ( 18
)
Where ∆t is a fixed simulation sampling time. The cartesian position in each step x k is then calculated by evaluating the robot's forward kinematics

x k = f DK (q k ). (19)
Each simulation roll-out considers a fixed joint torque vector τ and simulates the robot dynamics through the horizon time t h performing N = t h /∆t steps. Each vertex v i of reachable space polytope P x , calculated for a joint position q k and a horizon t h , is generated by a joint torque vector τ i , belonging to a subset of vertices of the joint torque polytope P τ .

To evaluate the accuracy of the proposed method, the difference in reached space produced by the constant liner model (polytope P x ) and time varying nonlinear model of the robot is evaluated, for the same set of applied joint torques. The applied joint torques chosen for the experiments are the joint torques τ i generating the vertices v i of the reachable space polytope P x calculated by the linear model.

As shown on figure 5, for each of the n v vertices v i of the polytope P x the simulation (18) is performed and all the cartesian positions x k of the robot in each of the N = t h /∆t sample times are retained for the further analysis.

X = {x 0,0 . . . x 0,N τ0 , x 1,0 . . . x 1,N τ1 , . . . , x nv,0 . . . x nv,N τn v } + = = =
Fig. 6. Depiction of the defined quantitative metrics for a simple 2 link planar robot.

First metric, m 1 , is defined as a ratio of number of the simulated cartesian robot positions x k ∈ X inside the polytope P x with respect to the overall number of simulated positions.

m 1 = |X ∩ P x | |X | = |X ∩ P x | N • n v (20)
The metric m 1 gives an insight on how well the polytope P x encapsulates the real reachable space determined with simulations. This metric might be used to asses the confidence in the polytope P x for safety applications, the closer this metric is to 1, the less chance that the robot can exit the space bounded by P x .

The second metric, m 2 , is defined as the ratio of the volumes of the polytope P x and the convex hull CH of the points

x k ∈ X ∩ P x , R 1 = CH(X ∩ P x )
representing a subset of the points reached by the robot in the simulations found inside of the polytope P x .

m 2 = V ol{R 1 } V ol{P x } (21)
The metric m 2 gives an insight on how much volume of the polytope P x is actually attainable by the robot, with respect to the simulations. This metric is designed to asses the quality of the approximation, if the score on this metric is low, this means that the large part of the P x is not actually reachable. The third metric, m 3 , is defined as the ratio of the volumes of the polytope P x and the convex hull CH of all the points

x k ∈ X , R 2 = CH(X )
reached by the robot in the simulations.

m 3 = V ol{P x } V ol{R 2 } (22)
The metric m 3 gives insight about the ratio of the volumes of the polytope P x and the simulated reachable space of the robot CH(X ). This metric is designed to asses the quality of the volume of the P x as a metric. The further this metric from 1 the less confidence one has in the volume of the P x .

Benchmark comparison

To benchmark the accuracy of convex polytope approach P x , it is compared with the coarse approximation of the robot's end-effector reachable space using fixed cartesian velocity ẋ and acceleration ẍ limits specified by robot manufacturer.

ẍ ∈ [ ẍmin , ẍmax ], ẋ ∈ [ ẋmin , ẋmax ] (23)
The cartesian space limitations (23) assume that the robot's acceleration ẍ and velocity ẋ capacity is constant, which is a not true since the robot's capacity depends highly on robot's joint configuration q and velocity q [2] [START_REF] Skuric | On-line force capability evaluation based on efficient polytope vertex search[END_REF]. However, they are very convenient and widely used in many applications. A simple convex space C x (cube in 3d) of reachable end-effector positions given the cartesian space limitations (23) can be calculated

C x = { x k+1 ∈ R m | x k+1 = ẍk t 2 h 2 + ẋk t h + x k ẍk ∈ [ ẍmin , ẍmax ] ẍk t h ∈ [ ẋmin , ẋmax ] } (24)
Figure 7 shows the visual comparison of the two approaches: P x and C x calculated for one pose of the Franka Emika Panda robot.

Results

The proposed analysis of the reachable space polytope accuracy is illustrated on the simulation of the collaborative robotic manipulator Franka Emika Panda with 7 degrees of freedom. Robot modeling and kinematics as well as simulations of robot dynamics are built using the python implementation of the roboticstoolbox [START_REF] Corke | Not your grandmother's toolbox-the robotics toolbox reinvented for python[END_REF] package.

The evaluation of the reachable space polytope P x and C x is carried out using the ICHM algorithm's efficient python implementation within the pycapacity [START_REF] Skuric | A real-time capable task-space capacity calculation python package[END_REF] package. The ICHM's approximation accuracy used for all the evaluations is δ=1mm.

Eight different horizon lengths t h are chosen for the experiments t h = [0.05s, 0.15s, 0.25s, 0.5s, 0.75s, 1.0s, 1.5s, 2.0s] and the simulation sample time has been set to ∆t=5ms.

The proposed method implementation as well as the code used for the experiments is fully open-source and can be found in the gitlab repository 4 . The gitlab repository additionally contains a Robot operating system (ROS) implementation of the proposed method for more interactive evaluation.

All the simulations are run on a computer with 1.90GHz Intel i7-8650U processor and 32Gb of RAM memory.

Execution time

Figure 8 shows the polytope P x execution time averaged over 1000 random robot configurations for each of the 8 horizon lengths. From the figure it can be seen that the execution time is relatively consistent through all the horizon lengths, with the average value of around 50ms. The constant execution time is expected because the horizon length t h does not have a big influence on the geometrical complexity of the polytope (number of vertices and faces) but its scale, and the scale itself does have a big influence the execution time.

Figure 9 shows the ICHM execution time for the 8 horizon lengths with different numbers of environmental constraints A e x ≤ b e . All the environmental constraints are generated randomly and all the results are averaged over 1000 random robot configurations for each horizon length t h and all 4 different numbers of environmental constraints: 0 (no environment constraints), 10, 100, 500 and 1000. Results show that up to 10 added environmental constraints the execution time of the method does not change significantly, having the average execution time around 50ms. With 100 environmental constraints the average execution time increases around 20% to around 70ms, for 500 the average execution is 100ms and with 1000 constraints the average execution time more than doubles, to 150ms.

The execution time analysis of the ICHM algorithm has shown that the average polytope P x enumeration takes between 50 and 150ms, which opens many doors for potential real-time applications.

The average execution time of the nonlinear simulation for t h =50ms and ∆t=5ms is around 60s, and it grows to about 240s for t h =2s. The cube C x average execution time stays constant for all horizon lengths at around 2ms.

Accuracy analysis

Figure 10 shows the evolution of the three metrics m 1 , m 2 and m 3 with respect to the horizon length t h for the P x and C x , averaged over 100 random robot configurations.

The first graph on the figure 10 shows the metric m 1 . It can be seen that the ratio of simulated robot positions x k ∈ X inside the polytope P x decreases with the increasing horizon length t h , while it increases for the C x . For larger horizon lengths the robot is able to move far from the initial joint configuration q at the beginning of the horizon, making the assumption of linearity of the robot model, made for calculating P x , largely inaccurate. The approach C x , on the other hand, improves with the increase in horizon length, since the space bounded by the cube C x increases linearly with the horizon length t h , at some point it encapsulates the whole workspace of the robot. However, it can be seen that even though the ratio m 1 decreases, the polytope P x still contains the most (≥60%) of the simulated robot's reached positions x k ∈ X . Finally, it can also be seen that the variance of the m 1 increases considerably, lowering the confidence of the polytope P x for longer horizons.

Evolution of the metric m 2 is shown on the middle graph of the figure 10 . It can be seen that the volume of all the reached robot positions inside polytope P x and C x decreases considerably with the increase of the horizon time t h . For the horizon times under t h ≤ 0.25s the majority (≥50%) of the reachable space polytope P x has been reached by the robot simulations. However for the longer horizon times t h ≥ 0.5s the reached volume ratio m 2 drops under 20% of the total volume of the P x . In the case of C x , the majority (≥50%) of the volume of the C x is never reached by the robot for all tested horizon lengths t h .

The right graph of the figure 10 shows the evolution of the metric m 3 . The graph shows that the volume of the polytope P x and the cube C x becomes several times higher than the volume of the R 2 with the increase in the horizon time t h . For the horizon time of t h =2s the average volume of P x is 12 times larger than the volume of R 2 and the volume C x is more than 50 times larger than R 2 . At the same time the average reached volume R 1 inside the P x and C x consists of under 10% of their volume (metric m 2 ), which potentially makes both polytope P x and C x impractical for many of applications.

However, for the shorter horizon times t h ≤ 0.25s the ratio m 3 is very close to 1 for P x , which means that the reachable space polytope's P x volume corresponds well to the actual reached volume of the robot R 2 .

Overall, this numerical analysis shows that the reachable space approximation based on the convex polytope P x works reasonably well for low horizon times t h ≤ 250ms, while for the longer horizon times t h > 250ms the assumption of the linearity of the system produces a large error in the reachable space estimation, producing the polytope P x with volume several times higher than the real reachable space of the robot (shown by the metric m 3 ), most of which is not reachable by the robot at all (showed by the metric m 2 ). Finally, the results show that the polytope P x , even though a coarse approximation, has better performance across all the defined metrics then the C x based on the constant cartesian limits, at the same time having smaller average volume for the all the horizon times. These results remind, once again, that robot's cartesian capacity is not constant and by modeling their evolution, even in a coarse manner, a substantial improvement in the estimation accuracy can be reached.

Discussion

Approximating the reachable space of the robot using convex polytopes P x is an efficient and relatively fast way of gaining information about the potential robot's position for the horizon of interest. Even though this approximation approach is neither an under or an over approximation of the true reachable space of the robot, as shown in the numerical analysis, for shorter horizon times t h ≤ 250ms this metric is reasonably precise. Due to its efficiency, easy and intuitive integration of different environmental constraints, possibility to consider the robot's link geometry and the fact that its output has the form of a convex set of inequality constraints (convex polytope P x ), it has a potential to be a useful tool for robot performance and safety analysis. However there are several important limitations of this approach.

Main limitations

First of all, the reachable space of a robot is highly non-convex and non-linear. As shown in the numerical analysis section 4, approximating this space with convex polytopes is relatively precise only for small horizon times t h ≤ 250ms.

Secondly, robot dynamics is highly non-linear and time variant as well. When calculating the reachable space polytope P x , as formulated in section 2, the robot model is linearised in the robot's joint configuration q k and velocity qk at the beginning of the horizon and this linearised model is considered constant during the horizon length. The consideration of constant robot model and its linearity, is only valid for small robot displacements around the current robot state, which is again a strong assumption and only valid for short horizon times.

Finally, the polytope formulation from the section 2, considers only constant joint torques τ applied to the robot during the whole horizon length t h . The polytope P x calculates the joint torque vectors τ that will produce the largest possible displacements in cartesian space given the robot and environment constraints. This calculation is based on the robot's linearised model at the beginning of the horizon, which in term, once again, limits the precision of the approximation to the small displacements around the initial robot state, preventing longer horizon lengths t h , as confirmed in the analysis section 4.

Potential applications

One promising field of applications for the polytope based reachable space metric is the human-robot interaction domain. The polytope is essentially a triangulated mesh and can be easily visualised, which could be an easy and intuitive way to give insight to the human operators about robot's real-time capabilities. It has a potential to be used for the haptic control applications, as human operators might not be always aware of how close the robot is to the limitation of its capabilities or what is the possible space the robot can be in while teleoperating. Furthermore, this metric could be potentially used to to increase human operators awareness about the robot's state during the co-manipulation by visualising the polytope P x in real-time and even potentially increase the interaction safety.

The numerical analysis has shown that for shorter horizon times t h ≤ 250ms the volume of the polytope P x is comparable to the volume of the simulated reached space R 2 , as shown on the figure 10. The volume of the P x could potentially be used as a performance metric for robotic workspace analysis and for applications such as robot design.

Finally, as the polytope P x can be represented as a set of convex inequalities Hx ≤ d, this metric could potentially be used as a set of constraints for robot control. For example in case of the Model Predictive Control (MPC) where, in many cases, robot's dynamics is already considered linear and constant during given horizon time, which corresponds well with the assumptions of the polytope P x formulation.

Conclusion

In this paper a convex polytope based robot's reachable space approximation metric is introduced. The metric predicts the robot's reachable space for a given time horizon based on the robot's actuation capabilities and its kinematic constraints. The proposed metric can intuitively integrate the environmental constraints as well as the robot's link geometry. With the execution time of around 50ms, for standard 7dof collaborative robot, it has a potential to be used for real-time applications.

However, the metric relies on a linearised robot model, which limits its use to shorter horizon times. The analysis has shown that the accuracy of the approximation decreases significantly for horizon times longer than 250ms.

Fig. 1 .

 1 Fig.1. Three images on the left show the comparison of the size of the reachable space polytope Px of a Franka Emika Panda 3 robot's end-effector for 3 horizons (left to right) t h = 0.05, 0.15 and 0.25s, for the same configuration. Whereas the three figures on the right show the constraining effect of the joint position limits on the Px for a fixed horizon time (t h =0.15s). Third image from the right shows the robot in its initial position q = [0, 0, 0, -π/2, 0, 3π/5, 0], close to the center of the joint ranges. In the second image form the right, the robot's joints q1 and q3 at their limits, q = [0, -1.59, 0, -2.9, 0, 3π/5, 0, 0]. On the last image, robot's joints q0 and q2 are at their limits q = [-2.72, 0, -2.72, -π/2, 0, 13π/5, 0, 0], preventing the robot to rotate around the z axis in one direction.

Fig. 3 .

 3 Fig.3. Resulting reachable space polytope Px for a Franka Emika Panda robot's endeffector when integrating the environmental constraints, horizon time used is t h =0.15s. Environment is defined as z ≥ 0.5m and y ≥ -0.2m.

Fig. 4 .

 4 Fig. 4. First four images show the construction of the reachable space polytope of the panda robot's link 3, P xl3 . From left to right they show, robot's link 3 as a line segment (in red) with polytopes Pxs and Pxe, their convex-hull, an example of the enveloping space L3 of the link 3 (in red), and the final polytope P xl3 calculated as a convex hull of the reachable space polytopes of each one of the 8 vertices (shown in yellow) of L3. Last two images further show an example of 4 enveloping spaces Li for each one of the robot's links, and their reachable space polytopes P xli . Robot is in its initial configuration, and the horizon time used is 150ms.

Fig. 5 .

 5 Fig. 5. Images show the steps (left to right) of the nonlinear simulation procedure on a simplified 2 link planar robot. First the polytope Px of the end-effector is determined.For each vertex of the Px, joint torque τ generating the vertex is applied to the robot simulation (figures 2 and 3). All the robot's end-effector positions in all the simulation time steps ∆t are saved for further analysis.

Fig. 7 .

 7 Fig. 7. Comparison of the end-effector reachable space approximation using the convex polytope Px (blue) and the cube Cx (orange) for the Franka Emika Panda robot with the horizon time t h =200ms. They are both evaluated in robot's initial position and the figures show two different views.

Fig. 8 .

 8 Fig. 8. Execution time of the polytope Px enumeration using ICHM algorithm. Average and standard deviation calculated over 1000 executions for each of the 8 horizon lengths.

Fig. 9 .

 9 Fig. 9. Average execution time and standard deviation calculated over 1000 executions for each of the 8 horizon lengths, with different number of added environmental constraints, ranging from 0 to 1000.

Fig. 10 .

 10 Fig. 10. Accuracy measures calculated for Px and Cx, averaged over 100 random poses of the Franka Emika panda robot for each of 8 horizon lengths.

There is a coupling between the velocity and torque limits of an actuator related to its power. Yet, it is common practice to consider a subset of these limits such as torque and velocity can be chosen independantly one from the other.

More information about the Panda robot at https://frankaemika.github.io/docs/

https://gitlab.inria.fr/auctus-team/people/antunskuric/reachable_space/

This work was funded by the BPI France Lichie project.