Fenceless obstacle avoidance method for efficient and safe human–robot collaboration in a shared work space

In a given manufacturing setting where workers or robots are coexisting in a confined area and their movements are not coordinated due to loss in communication or because they are freely ranging relative to each other, the development of an onboard safeguard system for a robot becomes a necessity to reduce accidents while the production efficiency is uncompromised. This article develops a two-dimensional dynamics model that predicts the relative position between a robot’s end-of-arm tooling and an approaching object or threat. The safety strategy applied to the robot is derived from the calculation of three parameters: the time of collision predicted from the linear motion between the approaching object and the robot’s end-of-arm tooling, the relative absolute distance, and the overlapping area ratio. These parameters combined are updated in a cost function that is sufficiently alarming the collision severity of an approaching object in real time. This model enables deployment a safe and a productive collaborative interaction in the manufacturing environment where workers and robots are seemingly moving in close proximity within an open workspace with less safeguard barriers.


Introduction
There is a need to develop predictive collision methods to improve worker-robot safety due to the increase installation of industrial and service robotics system in factories worldwide, where 1.4 million new industrial robots were expected by 2019 attributable to the use of autonomous guided systems. 1 Over the years, assessment methods have been deployed to improve the safety and reliability of robots serving several industry sectors including production, warehouses and logistics, autonomous transportation, and health-care assistive mobile robots. 2 Occupational accidents have been reporting worker injuries-those due to engineering and human errors 3 such as lack or absence of safeguard system to stop or slow down the operation of robot when a worker accidently enters robots work envelop, modifications to machinery program, or bypassing safeguards as well as inexperience of workers. The occupational machine-related accidents accounted for over 5000 fatalities in the United States during 2003-2010, of which 23% of fatalities were due to stationary machinery and 77% caused by mobile machinery. Manufacturing and construction were among top hazard industries accounting for 14% and 22% of the fatalities, respectively. 4 In the United States alone, the Bureau of Labor Statistics revealed 717 fatal work injuries occurred during the year 2013 as a result of direct contact with object. 5 Reports in Canada and in the Netherland revealed that most accidents are caused by operators entering hazardous zone of machinery during maintenance and the handling of production disturbances. 6,7 The worker safety is typically monitored by speed and separation from a robot which together used to detect collision in advance by particular reaction strategies and safeguards tools. 8 The tools include not-contact ultrasonic position sensors, contact capacitive or piezoelectric force and pressure sensors, door interlock switch, virtual barriers enabled by laser curtain, restricted work envelop that allows program robot movement within the desired workspace, pressure mat to detect the presence of operator, optical cameras, motion capture systems, and many others. Moreover, one common fenceless technology is laser scanners which have been traditionally used to safeguard a predefined number of zones around dangerous machinery. 9 Position-based collision mitigation are applied today in cyber-physical system where robot workspace is typically predefined into three safety zones to detect the human presence or moving machinery like AGV's. These technologies focus on safety first regardless of the efficiency of the manufacturing process. For example, wearable surface pressure sensor array on collaborative robot, such as MRK-Systeme KR SI, Fanuc CR-35iA, ABB YuMi, UR5, and KUKA LBR iiwa may lead to a reduced robot speed or entire motion stop at object detection. 10 The common practice is to enable all safety tools regardless of its effect on the production rate. Of course, while priority number one is to protect people in workplace, it is not yet clearly understood the trade-offs between safety and efficiency of robot.

Related work
The production efficiency in advanced manufacturing environment might be improved through the increase of the number of collaborative robot installation, however, that comes at the cost of risking the safety of coexistence with the workers. Particularly, "just-in-time sensing" based safety methods might not be adequate to address the issue of downtime in production when accidents or interference have already occurred. On the other hand, predictive collision-free algorithms-such as dynamic programming, genetic algorithm, and distance transform methods 11 are commonly used to optimize the path and the speed of mobile robots in the presence of obstacles. Similarly, to improve the efficiency of the manufacturing and the safety between human and industrial robot arm in a shared workplace, the path of a worker roaming around end-of-arm tooling (EOAT) of robot can be predicted at an instance and kept updated in real time with goal to estimate the future level of collision before accidents happen.
The goal of robot path planning in dynamic environment is to create a collision-free route more accurately and timely. Masehian et Al. survived over 1381 articles covering classic and heuristic approaches in robot motion planning within time span 1973 to 2007. 12 Classic methods were popular prior to 1982 and it uses general approaches such as roadmap, cell decomposition, potential fields, and mathematical programming. Classical methods such as PF, DWA, A*, and PRM path planning are easy to implement, perform better when blended with other methods, fast with respect to response time. 13,14 However, such conventional methods are computationally expensive and typically rely on accuracy of cameras and sensors to collect data from the environment; making measurement uncertain when noise are inherent. 15 For example, a heuristic A-Star and dynamic steering algorithms use predefined grid-based environment map to find the shortest path based on information received from sonar sensor. 16 Distance remote control algorithm was developed to predict and avoid collisions between vehicles by maintaining a given safe distance between the robot and the obstacle. 17 Potential fields were first introduced by Oussama Khatib in 1986, where an attractive potential pulls the robot toward the goal configuration, and a repulsive potential would push the robot away from the obstacles. 18 The mathematical programming approach formulates an optimization from a set of inequalities on the configuration parameters that minimize collision. The dimension of the configuration space implemented in cell-decomposition is large and time-consuming. 19 Nikolakis et al. developed system that used depth data and tool center point coordinate system of a robot to calculate the distance between the robot's end-effectors and worker in the workspace to control the speed of the robot and send notifications or alerts. 20 Long et al. developed a three-dimensional (3D) dynamic security zone for a collaborative robot UR10 to increases its performance in terms of average velocity by using three post-collision safety laser system. 21 Object detection around this robot was divided into three security modes: nominal mode, reduced mode, and passive mode, which change automatically according to the distance. Chen et al. studied the interaction efficiency of mobile robot in human traffic area to guide future design of robot-assisted pedestrian evacuation algorithms. 22 In their research, a laser range finder was integrated in Adept Pioneer P3-DX mobile robot, and a social force model experiment was also applied. The results showed that robot can detect pedestrian and avoiding the collision while maintaining desired average speed. Other obstacle avoidance methods used in path planning of mobile robotics include obstacle-dependent Gaussian potential field. 23 This method is applicable in static and dynamic environment without consideration of acceleration or speed of mobile robot. Distance histogram Bug algorithm is based on average path length and mainly implemented for static environment. 24 Other emerging techniques such as neural network, Fuzzy logic, and genetic algorithms are used for complex situations such as uncertain workspace and multiple inputs/outputs. 25,26 Safety methods for human-robot interaction in a fenceless robotic cell have been investigated with effort to maximize admissible velocity of a robot and to reduce injury. The projection of human arm motion into robot path was investigated to estimate the collision of robot. 27 This method is based on the optimal path related to the shortest duration of movement using segmentation of path and nonlinear programming; however, it does not consider size of objects and how speed control can be implemented to improve safety. Moreover, it uses noncost effective motion-tracking system that might be impractical for industrial application. Polverini et al. introduced assessment method using Kinetostatic safety field index only to study real-time collision avoidance in unstructured and dynamic environments. 28 The method uses mesh technique to calculate the index based on relative position and velocity vectors in real time between a moving rigid body "source of danger" and a point in space. Ragaflia et al. introduced trajectory generation algorithm acquired from different depth sensors for estimation of the kinematic configuration of human worker moving inside the robotic cell to maximize productivity function subjected to safety constraint. 29 The proposed algorithm takes into account volumes as obstacles modify the path of preprogrammed task of robotic arm. Sharkawy et al. utilized torque sensors data in multilayer feedforward neural network to train the dynamic model of joint manipulator with and without external contacts. 30 Mainprice and Berenson presented a framework that allows the human and the robot to perform simultaneous manipulation safely based on early prediction of the human motion to enable collaborative tasks in close proximity. 31 Wang et al. introduced a safety index with a distance factor is defined to describe the level of safety during physical human-robot interaction. The optimization problem was solved within the sampling time of the system allowing algorithm be applied in real-time trajectory generation. 32,33 Kanazawa et al. obtained robot's trajectory based on a probabilistic prediction of the worker's motion and the reduction of waste time simultaneously. 34 Notable summary of advantage and disadvantage of popular obstacle avoidance method are presented in Table 1. More recent in-depth reviews of techniques related to the navigation of unmanned mobile vehicles through various environments with obstacles can be found in the following literatures. 35,36,37,38,39,40,41 This article presents the development of obstacle avoidance method for industrial robots in manufacturing workplace based on prediction of dynamic motion of an approaching worker or robot within a shared workplace. Furthermore, the article discusses a mechanism to maximize efficiency and safety of shared workspace. The article introduces real-time prediction model of collision between a stationary industrial robot and a moving target based on a cost function that accounts for estimated time of impact, distance to impact, and the affected area of collision. The article offers a novel continuous speed control of robot based on simple cost function. Unlike majority of traditional industry-based safety methods which relies on the events to take place or on enforcing robot to change its path, the goal is to warrant robot or worker enough time by slowing down robot's speed; therefore, giving worker enough time to take corrective action and avoid collision that may result in production downtime and injuries. The article is organized as follows: a description of the robotworker shared space is defined in the "Related work" section. The relative position derivation is obtained from Table 1. Summary of few path planning and obstacle avoidance methods.

Method
Advantage Disadvantage Probabilistic roadmaps High-speed, parallel implementation. 42 Path search for every individual query. 19

Simulated annealing
Deals with highly nonlinear models, chaotic and noisy data, and many constraints and can generate smooth path. 43 Very slow when cost function is expensive, cannot tell whether it has found an optimal solution. 44 Neural networks Suitable for artificial intelligence, can be trained to solve nonlinear static and dynamic without knowledge of environment structure. 45 Computationally efficient but depend on initial value. 46 In 3D environments, the computation complexity will explode simultaneously. 47

Genetic algorithms
Better optimality robust based on stochastic search. 48 Poor local optimizing, oscillation, and convergence speed largely depends on the encoding method. 49,14 Ant colony optimization Proved to be able to deal with multiobjective path planning problems, and it is also able to tackle continuous planning problem. 50 The efficiency yet needs to be tested in real physical experiment. 51 Particle swarm optimization Simple to implement and requires less computing time and performs well in various optimization problems. 52 Easy to fall into local optimum in high-dimensional space and has a low convergence rate in the iterative process. 53

Fuzzy logic
Ability to make inferences even in uncertain scenarios with static and dynamic obstacles. 54,55 Robot stops to update the environmental map to determine the next move. 56 distributed sensor measurement. Linear dynamic path prediction of worker is derived and used in collision and safety models. An illustrative simulation example is presented in detail in the "Method" section. Worker trajectory and robot path are used together to calculate collision parameters, which are then used in cost function. Finally, collision and safety model is suggested to calculate the probability of collision and its effect on robot efficiency.

Problem definition
This project investigates a work environment comprised of an industrial robot arm and a fast-moving worker shown in Figure 1. The goal is to track the position vector of the worker relative to the robot's EOAT and then predict whether the motion trajectories would lead to possible collision. A collision cost function that decides the threat level is defined here from three parameters, simultaneously: time-to-impact t c , the relative distance between the EOAT and worker D c , and the predicted overlapping area ratio at impact S c . The t c is calculated from the time it takes the linear trajectories of both EOAT and the worker to intersect. Therefore, the values of t c if exist vary from zero to infinity, where collision is more probable to occur if the t c value is close to zero. Likewise, a large value of t c detected within a small shared space refers to an insignificant relative speed measured between EOAT and worker. A false alarm could be generated in occasions when the value of D c is small but the relative speed is very small. In such scenario, it is important to introduce other indicators such as magnitude of overlapping between the two trajectories calculated at different times. Because the position vector of worker could be calculated in reference to the EOAT by using geometrical transformation, it becomes possible to define a dynamic safety area around robot's EOAT, say a circle or a square with characteristic length of R R , within which interference with worker would be restricted. Therefore, this allows selecting an area that encloses dangerous part of robot arm. Similarly, worker area is defined by characteristic length R T , and work envelope ⍀ e of a robot is defined by characteristic length R e . In later section, an energybased cost function is utilized to combine the aforementioned parameters. The value of such cost function could be evaluated in real time to evaluate the threat such as no intersection and worker is far-off robot's EOAT, that is, t c ¼ : , or a worker is marginally within work envelop of robot but the relative speed is small, that is, t c ) 0; R e > D c > R R ð Þ , or speed of worker is fast relative to robot and marginally close t c ! 0; R e > D c > R R ð Þ , or no intersection and worker is marginally close, that is, ðt c ¼ : f g; R e > D c > R R Þ, and many other combinations. Under this definition, the threat level is customizable to the safety and production requirements. Moreover, the location of the restricted area is dynamic and became a subset of the working envelops of the robot. Unlike majority of traditional methods which have deployed static zones around robot all times no matter where the EOAT is located, or in most ongoing research proposal for a zone with dynamic radius, our proposed method promote to utilize space effectively by alarming system around collision areas. Ultimately, this facilitates the installation large number of robots in a workplace.

Measurement of the relative position
A simplified two-dimensional (2D) kinematic model for the work environment in Figure 1 is obtained and shown in Figure 2. The dynamical model assumes two moving areas in planner workspace with the goal to better understand the trade-offs between safety and efficiency, while a 3D space model can be simply obtained by adding the third  independent motion into the model in terms of depth coordinate. An arrangement of at least two position sensors S 1 and S 2 would provide 2D measurements within a confined area. The sensors are fixed relative to a world Cartesian coordinate at absolute position vectors P 1o and P 2o , respectively. Network of sensors (such as laser, ultrasonic, and cameras) could be distributed and fused to measure larger area. In this research, the derivation presented for pair of sensors. A bigger area can be covered by adding replicas of this model all around the robot where each pair of sensor covering a segment. The absolute position vector of the EOAT P R o is recorded by the robot's encoders in real time t. The absolute position vector of the worker P T o is indirectly obtained by two linear distance measurements D T 1 and D T 2 . The position vector of the worker relative to the EOAT is P T R and equals to P T o À P R o . The derivation of the position vector of the worker follows the trigonometry in Figure 2. The angle formed by vectors P 1o and P 2o is simply given by Apply the cosine law for each of the two triangles formed by the origin of the world Cartesian coordinate system, the worker and each sensor. This yields The unknown distance P T o k k and the angle q 2 are obtained numerically by solving the two nonlinear equations simultaneously.
Then calculate the worker angle q T measured CCW from X-axis by using Finally, the position components of the worker relative to the world Cartesian coordinate are obtained as follows In summary, it is plausible to calculate the 2D position of the worker by using two noncollinear distance sensors. The accuracy of the worker's position is sensitive to the quality of the data measured in real time. Therefore, highend sensors (i.e. capable of extracting depth such as ultrasonic, laser, radar, optical camera, and motion capture system) are recommended together with classification and filtration techniques that help reduce undesired noises in the prediction model. For example, generalized hough transform technique could be used in vision systems to extract the key features of moving object including its dimension and distance. 57

Linear dynamic trajectory
Consider the worker position P T o is measured for t o À nDt t t o , where n is the number of samples, Dt is the sampling time increment of the sensors, and t o represents the present time or the initial time at which the future path of the worker is unknown. On the other hand, the industrial robot is preprogrammed and its EOAT position P R o is often registered and given for 8t ! t o . Because the worker position is unknown to the robot for t > t o , future collisions between the EOAT and worker cannot be identified unless there exists a trajectory model that helps predicts the future position of the workerP T o . The prediction trajectories create positions in area ⍀ p , where In this article, it is assumed that the dynamic motion of the worker at t ¼ t o is under constant force with constant acceleration vector € P T o t o ð Þ, and initial velocity vector This assumption could also be applied to the robot in the absence of encoder readings, where the predicted position of EOATP R o replaces the true value P Ro . Therefore, the predicted position vector of the worker can be written by using equation of motion at constant acceleratioñ where the prediction time t p ¼ 0; Dt f p Â Ã is simulated until a final time Dt f p with a desired sampling rate Dt p . The velocity and acceleration are calculated by using backward finite difference method 58 Therefore, it is sufficient to store three position samples at t ¼ t o ; t o À Dt ; t o À 2Dt f g to predict the trajectory of the worker. While the rate of change of worker's velocity is assumed to be constant at instant, the model keeps updating the value of acceleration regardless of its magnitude. It shall be noted here that the proposed trajectory model is also deployable in case of robot-robot interaction and might produce more accurate results specially because human motion could be based on random walk.
It should be noted that the trajectory in equation (6) provides a realistic assumption about the motion of the worker. This could be justified because the trajectory is periodically updated. Other models relevant to social behavior could be investigated such as random walk, Markov Chain, and swarm or Brownian motion. Moreover, the linear trajectory model should be sufficient to predict the motion of worker when sampling rate is increased; however, the risk of collision is expected to be higher.

Collision model
A collision occurs in an areas ⍀ c at time set t c ¼ t 1 ; . . . ; t k g f and when the predicted path of the workerP T o t ð Þ and the true path of the EOAT P R o t ð Þ intersect inside the workplace envelope, that is, ⍀ c & ⍀ e . The worker and the EOAT are defined here by 2-D areas such as a square or circular shape. For example, the parametric functions defining the position of circular worker area ⍀ T are computed for 8v and 8r T where v 2 0; 2p ½ and r T 2 0; R T ½ . Similarly, the parametric functions defining the position of the circular moving area ⍀ R of the EOAT are computed for 8v and 8r R where r R 2 0; R R ½ . Regardless of the shape of the areas representing the robot or worker, the collision occurs when the following criterion holds where Dt r f is the final real-time simulation used to study the robot and worker movement. The solution of the criteria provides a time set for collisions t c with corresponding locations. Therefore, the time-to-impact is defined here by the time it takes a worker or a robot to respond for such threat is t c À t o . In addition, the overlapping area ratio S c at impact t c is a useful measure of safety and it is defined here by the ratio of the overlapping areas to the total areas, that is, All of these predictive measures with the true distance between the centers of the worker and EOAT at t o , that is, D c ¼ P TR t o ð Þ k kare independent and can be used to develop strategies for mitigation the risk of collision. It should be noted here that those control variables comprise the collision model M c and must be all updated in real time as the position of the worker changes. Also, the computational efficiency required to solving equation (11) can be substantially improved by narrowing down the search window within the area range covered by sensors.

Safety model
The goal of a safeguard system is to continuously control the velocity _ P R o t i ð Þ of the robot at any sampled time t i , depending on the predicted level of collision. The severity of the collision is obtained from the safety control model M S t i ð Þ in Figure 3, which uses collision model variables as control inputs: (1) inverse of time-to-impact The constant values e and d are small real numbers of order magnitude between fÀ1 to À4g used to bound the value of the cost function from reaching infinity when the relative distance or the time to collision are zero. The purpose of x 4 is to supervise or adjust the response to desired standard practices. Therefore, the worker is considered to be safe when all the state variables go to zero, that is, x 1 ! 0 when t c À t o ð Þ!1 , and x 3 ! 0 when S c ! 0, x 1 ! 0 when t c À t o ð Þ!1 . The overlapping area ratio-if exists-is bounded S c 2 0; 1 ½ . According to the definition of the overlapping area ratio, its maximum value S cmax can be achieved when the worker area ⍀ T is completely inside the robot area ⍀ R or vice versa. Because the areas are already known, their values can be calculated as follows: While there are many unaccounted variables that effect the safety, it is crucial to design a safety model that controls the speed of the robot using input vector X. Figure 3 suggests a safety model comprised of cost function C and SISO controller M A . It is desired to produce a bounded output Y that is not only proportional to the cost function but also can be used to easily override the speed of the robot. For example, consider safety model M s be comprised of multiplication of C ¼ X W X T and attenuator M A ¼ 1=A, where W is 4 Â 4 matrix and A is real number. The cost function is studied for a desired collision area ⍀ m & ⍀ e . The safety model outputs become bounded quantity 1 ! Y ! 0. The upper values of Y indicate that the robot can be operating close to the maximum speed, and the lowest values imply that the robot speed should be slowed down.

Real-time tracking and calculation of true impact
In a real-time prediction model, a discrete-step solver for simulation could be chosen for robot controller. Most modern robot controllers are running a cyclic control software, for example, 10 ms (or the better one 1 ms) for the position control loop. KUKA, for example, provides an RSI interface running at 1 ms and FANUC has a similar system called "dynamic path modification." Furthermore, vision system suffers from delay time during tracking of an object. Vision systems typically do not operate at 1 kHz rather most vision system run at 10-50 frames/s. In this simulation, long-range laser triangulation distance sensor is assumed, which typically has measurement rate much better than 1 kHz.
In this example, a robot moves in circular coordinate with parametric functions described by The position of the worker is measured indirectly by two linear distance sensors at sampling rate Dt ¼ p =20. The sensor measurements are simulated where sensor 1 reads D T 1 ¼ 5 ffiffi t p and sensor 2 reads The position components of the worker relative to the world coordinate are reversed engineered by applying equations (1)(2)(3)(4)(5). The distance ||P T o || and the angle q 2 are solved for every time step using the "fsolve" function in MATLAB software. 59 Then the worker P T o ð Þ x ; P T o ð Þ y positions are computed and plot with robot position Figure 4. The robot area ⍀ R and the worker area ⍀ T are moving in spatial domain and simulated by square shapes whose lengths are r R ¼ 1 and r T ¼ 0:4, respectively.
A true or "physical" impact between the worker and the robot is said to occur when the true moving areas physically overlap, that is, the impact time t c 6 ¼ : (9) and then substituting the latter values in equation (11). The case study is studied for period of time t ¼ p =2; p ½ at sampling rate Dt p ¼ Dt . The overlaping condition is checked at every time step, and the collision variables are recorded in Table 2. The results show that the true impact time took place at t c ¼ 2:1991 where the worker absolute position is P T o ¼ 9:6317. The percentage of the overlapped areas is S c ¼ 1:64% at a distance of P TR k k ¼ 0:4876 calculated between the center points of the worker and the robot. The results presented in Table 2 assume no prediction model enabled. The next section computes future collision from dynamic prediction equation of the worker.

Predictive impact model
Real-time tracking provides current and prior information for the position of the worker as function of time, but future positions are unknown due to the randomness of the human intention. Therefore, majority of safeguard technologies used in industrial robotics environment rely on the actual presence of the worker, where robot's controller is alerted only when the worker resides within predefined zones that are fixed around the space of the robot. Such alerts might be too late to prevent accidents especially when the worker or the robot movement is too fast! Instead, a time-varying predictive model of the worker path was introduced in equation (6) to provide a set of future positions of the worker calculated at every instantaneous position of the robot. This simply appears as trajectories tangent to the path of worker.
In the following simulation example, the case study in the previous section is used to predict the trajectory profile of the worker by using equation (6)(7)(8)(9) starting at absolute time instant t ¼ t o ¼ 6p =10¼1.885. The profile is generated for a period of time Dt f ¼ 2 and plot in Figure 5. Also, the true path of the robot is simultaneously simulated using given parametric functions. Criterion in equation (11) is evaluated at time t o and with worker trajectory generated for future time period t p ¼ 0; 2 ½ . Table 3 provides a summary of the simulation results and shows that the robot center is 1.1162 apart from worker center at time t ¼ t o ¼1.885 with estimated time-to-impact is 0.314. The percentage of area overlapped at impact moment is 2.4% with 0.4766 apart from each other's center points.  The above procedures are repeated for every time step t o in the time range t ¼ p =2; p ½ and with time step increment of Dt ¼ p=20. The first three increments are ignored because the prediction model requires prior knowledge of the past three positions. The results in Table 4 show prediction of two impacts estimated at times t o ¼ 1:885; 2:042 f g and with time-to-impacts are f0.314, 0.157g, respectively. The third predicted collision coincides with true impact where t o ¼ 2:199 f gand timeto-impact is zero.

Cost function and the probability of collision
Let the multi-input-single-output model M S ¼ X W X T =A be bounded in a spatial domain defined by the positions of a worker traveling along predicted paths shown in Figure 6 and with overall area defined by ⍀ m ¼ ⍀ p . The design of the coefficients in the weight matrix of the cost function should be such that it produces maximum values when the true collision condition is satisfied. The cost function should also provide an indication of collision severity proportional to the value of the maximum values. According to the simulation results in Table 4, there are three possible future collisions, with one being a real impact at t c À t o ð Þ%0. We want to inspect a diagonal weight matrix for the cost function C. This eliminates cross-correlation between the state variables and assumes orthogonality, that is, independent from each other. The following cost function is suggested where v is the coefficient corresponding to the state variable x. This cost function is computed for 8t o . In this special case study, it is assumed that there is no external control signal, that is, v 4 ¼ 0. The values of x 1 and x 3 are zero when there is no predicted collision. However, 0 < x 3 < S c max when a collision is predicted. In this simulation example, we where a is a correction gain used to amplify the effect of the overlapping term in the cost function. Let e ¼ Dt p =2 ¼ 0:0785, d ¼ 0:001 and a ¼ 100. The cost function is evaluated for a range of time t ¼ p =2; p ½ with step size Dt ¼ p=20. In general, when the cost function value is C ) 0, a collision is expected, while it becomes safe when the value C ! 0 as shown in Figure 7. More particularly, the set of solutions in Table 4 generates large values and with the cost function of the true impact being the maximum. The cost function response in Figure 7 has Gaussian distribution with one peak at the time of true impact t o ¼ 2:199 f g. The positive value of the slope d C =d t within left side of the bell shape in Figure 7time range p =2; 2:199 ½ -indicates an increasing risk of collision as the time progresses. This is because the worker Figure 5. The simulation of the predicted trajectory of the worker starts at t o ¼ 6p=10 for time period of Dt f ¼ 2 and with sampling rate of Dt p ¼ p=20. The path of the robot is generated for simultaneously using given parametric functions. The first three sampled steps are not included in prediction. Table 3. Predicted impact between worker and robot studied at  Figure 6. Simulation of the predicted trajectories of the worker and the true path of robot generated over t ¼ p=2; p ½ . International Journal of Advanced Robotic Systems is approaching the true impact condition. However, the negative slope in right side of bell shape-time range 2:199; p ½ -means that the worker is moving away from collision and therefore become safer. Moreover, one could set a threshold value based on Table 4 and above which all the predicted and true collisions would occur.
The cost function data in Figure 7 could be fitted into a combined Gaussian curve of form where k is the number of Gaussian functions. The secondorder Gaussian, that is, k ¼ 2 can be obtained by using MATLAB statistical toolbox "fit" function with "guess2" type, which yields C t ð Þ ¼ 3:41e The above model provides accurate fitting of cost function at any time. To better understand the meaning of the distribution function, a first-order Gaussian distribution is generated using "guess1" function in MATLAB and then plotted together with the second-order and the raw data in Figure 7. The fitted first-order Gaussian function with This first-order Gaussian distribution reveals that the true impact happens at the mean value ¼ 2:177 with standard deviation s ¼ 0:1909. The curve follows the normal distribution where the cumulative distribution function could be obtained by integrating equation (15). The "normcdf" function in MATLAB software calculates the cumulative distribution function evaluated at time t. Figure 8 is a plot of the collision probability P t ð Þ evaluated over time range [À1; t . Therefore, the probability of a collision in time range t 1 ; t 2 ½ , that is, t 2 ! P t ð Þ ! t 1 is equal to P t 2 Þ À P t 1 Þ ð ð . Alternatively, the probability of collision can be obtained by subtracting areas under the curve of Figure 8. The figure indicates that the collision is more likely to occur as human approach time .
Higher-order Gaussian combinations are suggested for scenarios involving multiple true impacts which appear in the form of peaks in the cost function curve. For example, k ¼ 4 could provide a good fit for two true impacts. However, other forms of cost function, different from quadratic form proposed in equation (12), could produce non-Gaussian distribution.

Safety control system
The controller model selected in this example is a simple attenuator 1=A where the value A ¼ max: C is calculated for all sum of areas navigated by the worker and the robot. In the previous example, the total area, shown in Figure 6, has minimum attenuation of 1=A ¼ 0:24. The value of C t ð Þ=A normalizes the safety model output to a scale of 100%. This guarantees that the output of the controller model is proportional to the risk of collision and, therefore, can be directly used to control the speed of the robot v R o t ð Þ. One could program the robot's joint or axial speed In the example studied previously, the robot moves in circular coordinate with of the robot corresponding to the threat posed by the worker movement is plot in Figure 9. As expected, the speed of the robot winds down as the worker approaches dangerous conditions, comes to full stop at collision moment, and then gains speed quickly as the worker moves away. It shall be noted here that the original path of robot is kept unchanged and the collision avoidance is presented in the sense of zero contact force at contact between the colliding objects.

Robot productivity
Assume the robot productivity index is measured by the total energy cost function computed over a range of time When there is no worker around the robot, the robot speed remains unchanged. In the example discussed earlier, an uncontrolled robot operates at normal speed of v R o ¼ 3, where its total energy over time p =2; p ½ is computed by using equation (16). The productivity index is equal to 3p =4. However, productivity index of a robot under controlled speed can be computed from the area under the curve in Figure 9 or by numerically integrating the productivity function Prod ¼ 1:5 1 À C t ð Þ=0:24Þ ð , with C t ð Þ being the function obtained in equation (14). The value of the productivity index under a controlled speed-the worker is present and the controller is enabled-being 1.778. The overall reduction in the robot productivity as results of applying the proposed safety control system is only %24.5. Moreover, another alternative productivity index analogous to energy function could be obtained by raising the power of the robot's speed in equation (16) by the value of 2.

Discussion
The collaborative human-robot interaction introduced in this research acquires information from the pair of sensors to update the robot Euclidean position in real time. In an unknown manufacturing environment, where worker's motion trajectory is uncertain, an auxiliary microprocessor predicts future positions of the worker by assuming a worker moving in a constant acceleration. Because the proposed method uses minimum number of sensors, high sensitive and accurate distance sensors should be installed to improve the reliability of the calculated safety level.
The trajectory generated by the microprocessor is kept updated in real time. On the other hand, because the robot path is typically preprogrammed and deterministic, it is compared with the generated trajectories in real time with goal to search for possible future intersection(s). The simplicity of the mathematical model made it easy to impose multiobjective functions and control scheme. This prediction model was then combined with manufacturing efficiency through a contentiously updated cost function which is directed to control the speed of the robot. The key approach for predicting a collision is pertained by applying four factors that encapsulate the predicted time to collision, distance to collision, overlapping area at collision, and a manual remote control parameter determined by the user.
The simulated case study produced a cost function with multimodal Gaussian distribution. The mean values of the distribution are located at the peaks and correspond to high probability of collision. At such instants, the robot smoothly decelerates to a quasi-static stop and then accelerates to continue its normal operation specially when the threat is sufficiently far-off. The smaller the variance of the distribution is, the less time is given to the robot to react to the threat. However and in accordance to the proposed productivity index-which was defined by area under speed and time curve-the robot efficiency increases when the gradient of the cost function is steep.
The physical realization of the proposed prediction method becomes feasible to implement mainly because it requires few measurements as compared to classical and modern obstacle avoidance methods. These data are collected from a pair of sensors configured to scan a sectional area whose size depends on the range of sensing technology used. The method can be extended to cover larger area by summing data from other pairs of distributed sensors.

Future work
The preliminary setup in Figure 10 is being constructed to deploy the developed safety algorithm in an open workspace environment where a threat-like a worker-is prominent. The system is comprised of Mitsubishi RH-6FHD robot with CR-750 controller, and two low-cost longrange ultrasonic sensors connected to an external Arduino controller. The worker position and the cost function are being calculated in real time using the Arduino microprocessor, which then sends signal to the robot I/O. Currently, the Arduino controller is programmed to produce three digital outputs, and a robot program-written in MELFA RT-toolbox software-operates at three distinct speeds depending on the value received from I/O. However, this setup requires preprogramming and synchronizing with the motion of the EOAT, which is impractical for realtime experiment and verification. Therefore, we plan to continue with development by enabling Arduino and robot controllers to talk back to each other using analog signal instead of digital. More particularly communicate the position of the EOAT to PC and have it scanned by Arduino via USB port.

Conclusion
We developed a continuous speed control system to improve human-robot safety while robot productivity is not compromised. The predictive dynamic motion model for a human moving relative to a robot provided an accurate mean to estimate the future impact conditions which are used to construct a Gaussian cost function. This is used to calculate the probability of collision within a range of time and design a safety model to control speed of the robot in real time and based on the overall threat level. The case study showed that for an EOAT of a robot moving in deterministic motion under no safety control, an arbitrarily approaching worker poses a threat that would lead to a physical impact. However, when the safety control was deployed, the robot's speed responded proportional to the cost function and led to only 24% reduction in productivity during entire operation. This is very promising result specially because it does not lead to complete downtime in robot operation. The proposed algorithms can be used in application involving object locomotion such as mobile robots, driverless cars, and unmanned aerial vehicles.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.